File indexing completed on 2024-10-08 23:09:57
0001 #include "RecoJets/JetProducers/interface/PileupJetIdAlgo.h"
0002
0003 #include "DataFormats/JetReco/interface/PFJet.h"
0004 #include "DataFormats/JetReco/interface/Jet.h"
0005 #include "DataFormats/VertexReco/interface/Vertex.h"
0006 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0007 #include "DataFormats/Math/interface/deltaR.h"
0008 #include "FWCore/ParameterSet/interface/FileInPath.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "CommonTools/MVAUtils/interface/GBRForestTools.h"
0011
0012 #include "TMatrixDSym.h"
0013 #include "TMatrixDSymEigen.h"
0014
0015 #include <utility>
0016
0017
0018 const float large_val = std::numeric_limits<float>::max();
0019
0020
0021
0022 PileupJetIdAlgo::AlgoGBRForestsAndConstants::AlgoGBRForestsAndConstants(edm::ParameterSet const& ps, bool runMvas)
0023 : cutBased_(ps.getParameter<bool>("cutBased")),
0024 etaBinnedWeights_(false),
0025 runMvas_(runMvas),
0026 nEtaBins_(0),
0027 label_(ps.getParameter<std::string>("label")),
0028 mvacut_{},
0029 rmsCut_{},
0030 betaStarCut_{} {
0031 std::vector<edm::FileInPath> tmvaEtaWeights;
0032 std::vector<std::string> tmvaSpectators;
0033 int version;
0034
0035 if (!cutBased_) {
0036 etaBinnedWeights_ = ps.getParameter<bool>("etaBinnedWeights");
0037 if (etaBinnedWeights_) {
0038 const std::vector<edm::ParameterSet>& trainings = ps.getParameter<std::vector<edm::ParameterSet> >("trainings");
0039 nEtaBins_ = ps.getParameter<int>("nEtaBins");
0040 for (int v = 0; v < nEtaBins_; v++) {
0041 tmvaEtaWeights.push_back(trainings.at(v).getParameter<edm::FileInPath>("tmvaWeights"));
0042 jEtaMin_.push_back(trainings.at(v).getParameter<double>("jEtaMin"));
0043 jEtaMax_.push_back(trainings.at(v).getParameter<double>("jEtaMax"));
0044 }
0045 for (int v = 0; v < nEtaBins_; v++) {
0046 tmvaEtaVariables_.push_back(trainings.at(v).getParameter<std::vector<std::string> >("tmvaVariables"));
0047 }
0048 } else {
0049 tmvaVariables_ = ps.getParameter<std::vector<std::string> >("tmvaVariables");
0050 }
0051 tmvaMethod_ = ps.getParameter<std::string>("tmvaMethod");
0052 tmvaSpectators = ps.getParameter<std::vector<std::string> >("tmvaSpectators");
0053 version = ps.getParameter<int>("version");
0054 } else {
0055 version = USER;
0056 }
0057
0058 edm::ParameterSet jetConfig = ps.getParameter<edm::ParameterSet>("JetIdParams");
0059 for (int i0 = 0; i0 < 3; i0++) {
0060 std::string lCutType = "Tight";
0061 if (i0 == PileupJetIdentifier::kMedium)
0062 lCutType = "Medium";
0063 if (i0 == PileupJetIdentifier::kLoose)
0064 lCutType = "Loose";
0065 int nCut = 1;
0066 if (cutBased_)
0067 nCut++;
0068 for (int i1 = 0; i1 < nCut; i1++) {
0069 std::string lFullCutType = lCutType;
0070 if (cutBased_ && i1 == 0)
0071 lFullCutType = "BetaStar" + lCutType;
0072 if (cutBased_ && i1 == 1)
0073 lFullCutType = "RMS" + lCutType;
0074 std::vector<double> pt010 = jetConfig.getParameter<std::vector<double> >(("Pt010_" + lFullCutType).c_str());
0075 std::vector<double> pt1020 = jetConfig.getParameter<std::vector<double> >(("Pt1020_" + lFullCutType).c_str());
0076 std::vector<double> pt2030 = jetConfig.getParameter<std::vector<double> >(("Pt2030_" + lFullCutType).c_str());
0077 std::vector<double> pt3040 = jetConfig.getParameter<std::vector<double> >(("Pt3040_" + lFullCutType).c_str());
0078 std::vector<double> pt4050 = jetConfig.getParameter<std::vector<double> >(("Pt4050_" + lFullCutType).c_str());
0079 if (!cutBased_) {
0080 for (int i2 = 0; i2 < 4; i2++)
0081 mvacut_[i0][0][i2] = pt010[i2];
0082 for (int i2 = 0; i2 < 4; i2++)
0083 mvacut_[i0][1][i2] = pt1020[i2];
0084 for (int i2 = 0; i2 < 4; i2++)
0085 mvacut_[i0][2][i2] = pt2030[i2];
0086 for (int i2 = 0; i2 < 4; i2++)
0087 mvacut_[i0][3][i2] = pt3040[i2];
0088 for (int i2 = 0; i2 < 4; i2++)
0089 mvacut_[i0][4][i2] = pt4050[i2];
0090 }
0091 if (cutBased_ && i1 == 0) {
0092 for (int i2 = 0; i2 < 4; i2++)
0093 betaStarCut_[i0][0][i2] = pt010[i2];
0094 for (int i2 = 0; i2 < 4; i2++)
0095 betaStarCut_[i0][1][i2] = pt1020[i2];
0096 for (int i2 = 0; i2 < 4; i2++)
0097 betaStarCut_[i0][2][i2] = pt2030[i2];
0098 for (int i2 = 0; i2 < 4; i2++)
0099 betaStarCut_[i0][3][i2] = pt3040[i2];
0100 for (int i2 = 0; i2 < 4; i2++)
0101 betaStarCut_[i0][4][i2] = pt4050[i2];
0102 }
0103 if (cutBased_ && i1 == 1) {
0104 for (int i2 = 0; i2 < 4; i2++)
0105 rmsCut_[i0][0][i2] = pt010[i2];
0106 for (int i2 = 0; i2 < 4; i2++)
0107 rmsCut_[i0][1][i2] = pt1020[i2];
0108 for (int i2 = 0; i2 < 4; i2++)
0109 rmsCut_[i0][2][i2] = pt2030[i2];
0110 for (int i2 = 0; i2 < 4; i2++)
0111 rmsCut_[i0][3][i2] = pt3040[i2];
0112 for (int i2 = 0; i2 < 4; i2++)
0113 rmsCut_[i0][4][i2] = pt4050[i2];
0114 }
0115 }
0116 }
0117
0118 if (!cutBased_) {
0119 assert(tmvaMethod_.empty() || ((!tmvaVariables_.empty() || (!tmvaEtaVariables_.empty())) && version == USER));
0120 }
0121
0122 if ((!cutBased_) && (runMvas_)) {
0123 if (etaBinnedWeights_) {
0124 for (int v = 0; v < nEtaBins_; v++) {
0125 etaReader_.push_back(createGBRForest(tmvaEtaWeights.at(v)));
0126 }
0127 } else {
0128 reader_ = createGBRForest(ps.getParameter<std::string>("tmvaWeights"));
0129 }
0130 }
0131 }
0132
0133 PileupJetIdAlgo::PileupJetIdAlgo(AlgoGBRForestsAndConstants const* cache) : cache_(cache) { initVariables(); }
0134
0135
0136 PileupJetIdAlgo::~PileupJetIdAlgo() {}
0137
0138
0139 void assign(const std::vector<float>& vec, float& a, float& b, float& c, float& d) {
0140 size_t sz = vec.size();
0141 a = (sz > 0 ? vec[0] : 0.);
0142 b = (sz > 1 ? vec[1] : 0.);
0143 c = (sz > 2 ? vec[2] : 0.);
0144 d = (sz > 3 ? vec[3] : 0.);
0145 }
0146
0147 void setPtEtaPhi(const reco::Candidate& p, float& pt, float& eta, float& phi) {
0148 pt = p.pt();
0149 eta = p.eta();
0150 phi = p.phi();
0151 }
0152
0153
0154 void PileupJetIdAlgo::set(const PileupJetIdentifier& id) { internalId_ = id; }
0155
0156
0157
0158 float PileupJetIdAlgo::getMVAval(const std::vector<std::string>& varList,
0159 const std::unique_ptr<const GBRForest>& reader) {
0160 std::vector<float> vars;
0161 for (std::vector<std::string>::const_iterator it = varList.begin(); it != varList.end(); ++it) {
0162 std::pair<float*, float> var = variables_.at(*it);
0163 vars.push_back(*var.first);
0164 }
0165 return reader->GetClassifier(vars.data());
0166 }
0167
0168 void PileupJetIdAlgo::runMva() {
0169 if (cache_->cutBased()) {
0170 internalId_.idFlag_ = computeCutIDflag(
0171 internalId_.betaStarClassic_, internalId_.dR2Mean_, internalId_.nvtx_, internalId_.jetPt_, internalId_.jetEta_);
0172 } else {
0173 if (std::abs(internalId_.jetEta_) >= 5.0) {
0174 internalId_.mva_ = -2.;
0175 } else {
0176 if (cache_->etaBinnedWeights()) {
0177 if (std::abs(internalId_.jetEta_) > cache_->jEtaMax().at(cache_->nEtaBins() - 1)) {
0178 internalId_.mva_ = -2.;
0179 } else {
0180 for (int v = 0; v < cache_->nEtaBins(); v++) {
0181 if (std::abs(internalId_.jetEta_) >= cache_->jEtaMin().at(v) &&
0182 std::abs(internalId_.jetEta_) < cache_->jEtaMax().at(v)) {
0183 internalId_.mva_ = getMVAval(cache_->tmvaEtaVariables().at(v), cache_->etaReader().at(v));
0184 break;
0185 }
0186 }
0187 }
0188 } else {
0189 internalId_.mva_ = getMVAval(cache_->tmvaVariables(), cache_->reader());
0190 }
0191 }
0192 internalId_.idFlag_ = computeIDflag(internalId_.mva_, internalId_.jetPt_, internalId_.jetEta_);
0193 }
0194 }
0195
0196
0197 std::pair<int, int> PileupJetIdAlgo::getJetIdKey(float jetPt, float jetEta) {
0198 int ptId = 0;
0199 if (jetPt >= 10 && jetPt < 20)
0200 ptId = 1;
0201 if (jetPt >= 20 && jetPt < 30)
0202 ptId = 2;
0203 if (jetPt >= 30 && jetPt < 40)
0204 ptId = 3;
0205 if (jetPt >= 40)
0206 ptId = 4;
0207
0208 int etaId = 0;
0209 if (std::abs(jetEta) >= 2.5 && std::abs(jetEta) < 2.75)
0210 etaId = 1;
0211 if (std::abs(jetEta) >= 2.75 && std::abs(jetEta) < 3.0)
0212 etaId = 2;
0213 if (std::abs(jetEta) >= 3.0 && std::abs(jetEta) < 5.0)
0214 etaId = 3;
0215
0216 return std::pair<int, int>(ptId, etaId);
0217 }
0218
0219 int PileupJetIdAlgo::computeCutIDflag(float betaStarClassic, float dR2Mean, float nvtx, float jetPt, float jetEta) {
0220 std::pair<int, int> jetIdKey = getJetIdKey(jetPt, jetEta);
0221 float betaStarModified = betaStarClassic / log(nvtx - 0.64);
0222 int idFlag(0);
0223 if (betaStarModified < cache_->betaStarCut()[PileupJetIdentifier::kTight][jetIdKey.first][jetIdKey.second] &&
0224 dR2Mean < cache_->rmsCut()[PileupJetIdentifier::kTight][jetIdKey.first][jetIdKey.second])
0225 idFlag += 1 << PileupJetIdentifier::kTight;
0226
0227 if (betaStarModified < cache_->betaStarCut()[PileupJetIdentifier::kMedium][jetIdKey.first][jetIdKey.second] &&
0228 dR2Mean < cache_->rmsCut()[PileupJetIdentifier::kMedium][jetIdKey.first][jetIdKey.second])
0229 idFlag += 1 << PileupJetIdentifier::kMedium;
0230
0231 if (betaStarModified < cache_->betaStarCut()[PileupJetIdentifier::kLoose][jetIdKey.first][jetIdKey.second] &&
0232 dR2Mean < cache_->rmsCut()[PileupJetIdentifier::kLoose][jetIdKey.first][jetIdKey.second])
0233 idFlag += 1 << PileupJetIdentifier::kLoose;
0234 return idFlag;
0235 }
0236
0237 int PileupJetIdAlgo::computeIDflag(float mva, float jetPt, float jetEta) {
0238 std::pair<int, int> jetIdKey = getJetIdKey(jetPt, jetEta);
0239 return computeIDflag(mva, jetIdKey.first, jetIdKey.second);
0240 }
0241
0242
0243 int PileupJetIdAlgo::computeIDflag(float mva, int ptId, int etaId) {
0244 int idFlag(0);
0245 if (mva > cache_->mvacut()[PileupJetIdentifier::kTight][ptId][etaId])
0246 idFlag += 1 << PileupJetIdentifier::kTight;
0247 if (mva > cache_->mvacut()[PileupJetIdentifier::kMedium][ptId][etaId])
0248 idFlag += 1 << PileupJetIdentifier::kMedium;
0249 if (mva > cache_->mvacut()[PileupJetIdentifier::kLoose][ptId][etaId])
0250 idFlag += 1 << PileupJetIdentifier::kLoose;
0251 return idFlag;
0252 }
0253
0254
0255 PileupJetIdentifier PileupJetIdAlgo::computeMva() {
0256 runMva();
0257 return PileupJetIdentifier(internalId_);
0258 }
0259
0260
0261 PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet,
0262 float jec,
0263 const reco::Vertex* vtx,
0264 const reco::VertexCollection& allvtx,
0265 double rho,
0266 edm::ValueMap<float>& constituentWeights,
0267 bool applyConstituentWeight) {
0268
0269 resetVariables();
0270
0271
0272 const pat::Jet* patjet = dynamic_cast<const pat::Jet*>(jet);
0273 const reco::PFJet* pfjet = dynamic_cast<const reco::PFJet*>(jet);
0274 assert(patjet != nullptr || pfjet != nullptr);
0275 if (patjet != nullptr && jec == 0.) {
0276 jec = patjet->pt() / patjet->correctedJet(0).pt();
0277 }
0278 if (jec <= 0.) {
0279 jec = 1.;
0280 }
0281
0282 const reco::Candidate *lLead = nullptr, *lSecond = nullptr, *lLeadNeut = nullptr, *lLeadEm = nullptr,
0283 *lLeadCh = nullptr, *lTrail = nullptr;
0284
0285 std::vector<float> frac, fracCh, fracEm, fracNeut;
0286 float cones[] = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7};
0287 size_t ncones = sizeof(cones) / sizeof(float);
0288 float* coneFracs[] = {&internalId_.frac01_,
0289 &internalId_.frac02_,
0290 &internalId_.frac03_,
0291 &internalId_.frac04_,
0292 &internalId_.frac05_,
0293 &internalId_.frac06_,
0294 &internalId_.frac07_};
0295 float* coneEmFracs[] = {&internalId_.emFrac01_,
0296 &internalId_.emFrac02_,
0297 &internalId_.emFrac03_,
0298 &internalId_.emFrac04_,
0299 &internalId_.emFrac05_,
0300 &internalId_.emFrac06_,
0301 &internalId_.emFrac07_};
0302 float* coneNeutFracs[] = {&internalId_.neutFrac01_,
0303 &internalId_.neutFrac02_,
0304 &internalId_.neutFrac03_,
0305 &internalId_.neutFrac04_,
0306 &internalId_.neutFrac05_,
0307 &internalId_.neutFrac06_,
0308 &internalId_.neutFrac07_};
0309 float* coneChFracs[] = {&internalId_.chFrac01_,
0310 &internalId_.chFrac02_,
0311 &internalId_.chFrac03_,
0312 &internalId_.chFrac04_,
0313 &internalId_.chFrac05_,
0314 &internalId_.chFrac06_,
0315 &internalId_.chFrac07_};
0316 TMatrixDSym covMatrix(2);
0317 covMatrix = 0.;
0318 float jetPt = jet->pt() / jec;
0319 float sumPt = 0., sumPt2 = 0., sumTkPt = 0., sumPtCh = 0, sumPtNe = 0;
0320 float multNeut = 0.0;
0321 float sumW2(0.0);
0322 float sum_deta(0.0), sum_dphi(0.0);
0323 float ave_deta(0.0), ave_dphi(0.0);
0324 setPtEtaPhi(
0325 *jet, internalId_.jetPt_, internalId_.jetEta_, internalId_.jetPhi_);
0326 internalId_.jetM_ = jet->mass();
0327 internalId_.nvtx_ = allvtx.size();
0328 internalId_.rho_ = rho;
0329
0330 float dRmin(1000);
0331 float LeadcandWeight = 1.0;
0332 float SecondcandWeight = 1.0;
0333 float LeadNeutcandWeight = 1.0;
0334 float LeadEmcandWeight = 1.0;
0335 float LeadChcandWeight = 1.0;
0336 float TrailcandWeight = 1.0;
0337
0338 for (unsigned i = 0; i < jet->numberOfSourceCandidatePtrs(); ++i) {
0339 reco::CandidatePtr pfJetConstituent = jet->sourceCandidatePtr(i);
0340 const reco::Candidate* icand = pfJetConstituent.get();
0341 assert(icand);
0342 const pat::PackedCandidate* lPack = dynamic_cast<const pat::PackedCandidate*>(icand);
0343 const reco::PFCandidate* lPF = dynamic_cast<const reco::PFCandidate*>(icand);
0344 bool isPacked = true;
0345 if (lPack == nullptr) {
0346 isPacked = false;
0347 }
0348
0349 float candWeight = 1.0;
0350 if (applyConstituentWeight) {
0351 candWeight = constituentWeights[jet->sourceCandidatePtr(i)];
0352 }
0353 float candPt = (icand->pt()) * candWeight;
0354 float candPtFrac = candPt / jetPt;
0355 float candDr = reco::deltaR(*icand, *jet);
0356 float candDeta = icand->eta() - jet->eta();
0357 float candDphi = reco::deltaPhi(*icand, *jet);
0358 float candPtDr = candPt * candDr;
0359 size_t icone = std::lower_bound(&cones[0], &cones[ncones], candDr) - &cones[0];
0360
0361 if (candDr < dRmin)
0362 dRmin = candDr;
0363
0364
0365 if (lLead == nullptr || candPt > (lLead->pt()) * LeadcandWeight) {
0366 lSecond = lLead;
0367 SecondcandWeight = LeadcandWeight;
0368 lLead = icand;
0369 if (applyConstituentWeight) {
0370 LeadcandWeight = constituentWeights[jet->sourceCandidatePtr(i)];
0371 }
0372 } else if ((lSecond == nullptr || candPt > (lSecond->pt()) * SecondcandWeight) &&
0373 (candPt < (lLead->pt()) * LeadcandWeight)) {
0374 lSecond = icand;
0375 if (applyConstituentWeight) {
0376 SecondcandWeight = constituentWeights[jet->sourceCandidatePtr(i)];
0377 }
0378 }
0379
0380
0381 internalId_.dRMean_ += candPtDr;
0382 internalId_.dR2Mean_ += candPtDr * candPtDr;
0383 covMatrix(0, 0) += candPt * candPt * candDeta * candDeta;
0384 covMatrix(0, 1) += candPt * candPt * candDeta * candDphi;
0385 covMatrix(1, 1) += candPt * candPt * candDphi * candDphi;
0386 internalId_.ptD_ += candPt * candPt;
0387 sumPt += candPt;
0388 sumPt2 += candPt * candPt;
0389
0390
0391 frac.push_back(candPtFrac);
0392
0393 if (icone < ncones) {
0394 *coneFracs[icone] += candPt;
0395 }
0396
0397
0398 if (abs(icand->pdgId()) == 130) {
0399 if (lLeadNeut == nullptr || candPt > (lLeadNeut->pt()) * LeadNeutcandWeight) {
0400 lLeadNeut = icand;
0401 if (applyConstituentWeight) {
0402 LeadNeutcandWeight = constituentWeights[jet->sourceCandidatePtr(i)];
0403 }
0404 }
0405
0406 internalId_.dRMeanNeut_ += candPtDr;
0407 fracNeut.push_back(candPtFrac);
0408 if (icone < ncones) {
0409 *coneNeutFracs[icone] += candPt;
0410 }
0411 internalId_.ptDNe_ += candPt * candPt;
0412 sumPtNe += candPt;
0413 multNeut += candWeight;
0414 }
0415
0416
0417 if (icand->pdgId() == 22) {
0418 if (lLeadEm == nullptr || candPt > (lLeadEm->pt()) * LeadEmcandWeight) {
0419 lLeadEm = icand;
0420 if (applyConstituentWeight) {
0421 LeadEmcandWeight = constituentWeights[jet->sourceCandidatePtr(i)];
0422 }
0423 }
0424 internalId_.dRMeanEm_ += candPtDr;
0425 fracEm.push_back(candPtFrac);
0426 if (icone < ncones) {
0427 *coneEmFracs[icone] += candPt;
0428 }
0429 internalId_.ptDNe_ += candPt * candPt;
0430 sumPtNe += candPt;
0431 multNeut += candWeight;
0432 }
0433
0434 if ((abs(icand->pdgId()) == 1) || (abs(icand->pdgId()) == 2))
0435 multNeut += candWeight;
0436
0437
0438 if (icand->charge() != 0) {
0439 if (lLeadCh == nullptr || candPt > (lLeadCh->pt()) * LeadChcandWeight) {
0440 lLeadCh = icand;
0441 if (applyConstituentWeight) {
0442 LeadChcandWeight = constituentWeights[jet->sourceCandidatePtr(i)];
0443 }
0444 const reco::Track* pfTrk = icand->bestTrack();
0445 if (lPF && std::abs(icand->pdgId()) == 13 && pfTrk == nullptr) {
0446 reco::MuonRef lmuRef = lPF->muonRef();
0447 if (lmuRef.isNonnull()) {
0448 const reco::Muon& lmu = *lmuRef.get();
0449 pfTrk = lmu.bestTrack();
0450 edm::LogWarning("BadMuon")
0451 << "Found a PFCandidate muon without a trackRef: falling back to Muon::bestTrack ";
0452 }
0453 }
0454 if (pfTrk == nullptr) {
0455
0456 if (isPacked) {
0457 internalId_.d0_ = std::abs(lPack->dxy(vtx->position()));
0458 internalId_.dZ_ = std::abs(lPack->dz(vtx->position()));
0459 } else if (lPF != nullptr) {
0460 pfTrk = (lPF->trackRef().get() == nullptr) ? lPF->gsfTrackRef().get() : lPF->trackRef().get();
0461 internalId_.d0_ = std::abs(pfTrk->dxy(vtx->position()));
0462 internalId_.dZ_ = std::abs(pfTrk->dz(vtx->position()));
0463 }
0464 } else {
0465 internalId_.d0_ = std::abs(pfTrk->dxy(vtx->position()));
0466 internalId_.dZ_ = std::abs(pfTrk->dz(vtx->position()));
0467 }
0468 }
0469 internalId_.dRMeanCh_ += candPtDr;
0470 internalId_.ptDCh_ += candPt * candPt;
0471 fracCh.push_back(candPtFrac);
0472 if (icone < ncones) {
0473 *coneChFracs[icone] += candPt;
0474 }
0475 sumPtCh += candPt;
0476 }
0477
0478 if (icand->charge() != 0) {
0479 if (!isPacked) {
0480 if (lPF->trackRef().isNonnull()) {
0481 float tkpt = candPt;
0482 sumTkPt += tkpt;
0483
0484 bool inVtx0 = vtx->trackWeight(lPF->trackRef()) > 0;
0485
0486 bool inAnyOther = false;
0487
0488 double dZ0 = std::abs(lPF->trackRef()->dz(vtx->position()));
0489 double dZ = dZ0;
0490 for (reco::VertexCollection::const_iterator vi = allvtx.begin(); vi != allvtx.end(); ++vi) {
0491 const reco::Vertex& iv = *vi;
0492 if (iv.isFake() || iv.ndof() < 4) {
0493 continue;
0494 }
0495
0496 bool isVtx0 = (iv.position() - vtx->position()).r() < 0.02;
0497
0498
0499 if (!isVtx0 && !inAnyOther) {
0500 inAnyOther = vtx->trackWeight(lPF->trackRef()) <= 0;
0501 }
0502
0503 dZ = std::min(dZ, std::abs(lPF->trackRef()->dz(iv.position())));
0504 }
0505
0506 if (inVtx0 && !inAnyOther) {
0507 internalId_.betaClassic_ += tkpt;
0508 } else if (!inVtx0 && inAnyOther) {
0509 internalId_.betaStarClassic_ += tkpt;
0510 }
0511
0512 if (dZ0 < 0.2) {
0513 internalId_.beta_ += tkpt;
0514 } else if (dZ < 0.2) {
0515 internalId_.betaStar_ += tkpt;
0516 }
0517 }
0518 } else {
0519 float tkpt = candPt;
0520 sumTkPt += tkpt;
0521 bool inVtx0 = false;
0522 bool inVtxOther = false;
0523 double dZ0 = 9999.;
0524 double dZ_tmp = 9999.;
0525 for (unsigned vtx_i = 0; vtx_i < allvtx.size(); vtx_i++) {
0526 const auto& iv = allvtx[vtx_i];
0527
0528 if (iv.isFake())
0529 continue;
0530
0531
0532 bool isVtx0 = (iv.position() - vtx->position()).r() < 0.02;
0533
0534 if (isVtx0) {
0535 if (lPack->fromPV(vtx_i) == pat::PackedCandidate::PVUsedInFit)
0536 inVtx0 = true;
0537 if (lPack->fromPV(vtx_i) == 0)
0538 inVtxOther = true;
0539 dZ0 = lPack->dz(iv.position());
0540 }
0541
0542 if (fabs(lPack->dz(iv.position())) < fabs(dZ_tmp)) {
0543 dZ_tmp = lPack->dz(iv.position());
0544 }
0545 }
0546 if (inVtx0) {
0547 internalId_.betaClassic_ += tkpt;
0548 } else if (inVtxOther) {
0549 internalId_.betaStarClassic_ += tkpt;
0550 }
0551 if (fabs(dZ0) < 0.2) {
0552 internalId_.beta_ += tkpt;
0553 } else if (fabs(dZ_tmp) < 0.2) {
0554 internalId_.betaStar_ += tkpt;
0555 }
0556 }
0557 }
0558
0559
0560 if (lTrail == nullptr || candPt < (lTrail->pt()) * TrailcandWeight) {
0561 lTrail = icand;
0562 if (applyConstituentWeight) {
0563 TrailcandWeight = constituentWeights[jet->sourceCandidatePtr(i)];
0564 }
0565 }
0566
0567
0568
0569 float weight2 = candPt * candPt;
0570 sumW2 += weight2;
0571 float deta = icand->eta() - jet->eta();
0572 float dphi = reco::deltaPhi(*icand, *jet);
0573 sum_deta += deta * weight2;
0574 sum_dphi += dphi * weight2;
0575 }
0576 if (sumW2 > 0) {
0577 ave_deta = sum_deta / sumW2;
0578 ave_dphi = sum_dphi / sumW2;
0579 }
0580
0581
0582
0583
0584 assert(!(lLead == nullptr));
0585 internalId_.leadPt_ = lLead->pt() * LeadcandWeight;
0586 internalId_.leadEta_ = lLead->eta();
0587 internalId_.leadPhi_ = lLead->phi();
0588
0589 if (lSecond != nullptr) {
0590 internalId_.secondPt_ = lSecond->pt() * SecondcandWeight;
0591 internalId_.secondEta_ = lSecond->eta();
0592 internalId_.secondPhi_ = lSecond->phi();
0593 } else {
0594 internalId_.secondPt_ = 0.0;
0595 internalId_.secondEta_ = large_val;
0596 internalId_.secondPhi_ = large_val;
0597 }
0598
0599 if (lLeadNeut != nullptr) {
0600 internalId_.leadNeutPt_ = lLeadNeut->pt() * LeadNeutcandWeight;
0601 internalId_.leadNeutEta_ = lLeadNeut->eta();
0602 internalId_.leadNeutPhi_ = lLeadNeut->phi();
0603 } else {
0604 internalId_.leadNeutPt_ = 0.0;
0605 internalId_.leadNeutEta_ = large_val;
0606 internalId_.leadNeutPhi_ = large_val;
0607 }
0608
0609 if (lLeadEm != nullptr) {
0610 internalId_.leadEmPt_ = lLeadEm->pt() * LeadEmcandWeight;
0611 internalId_.leadEmEta_ = lLeadEm->eta();
0612 internalId_.leadEmPhi_ = lLeadEm->phi();
0613 } else {
0614 internalId_.leadEmPt_ = 0.0;
0615 internalId_.leadEmEta_ = large_val;
0616 internalId_.leadEmPhi_ = large_val;
0617 }
0618
0619 if (lLeadCh != nullptr) {
0620 internalId_.leadChPt_ = lLeadCh->pt() * LeadChcandWeight;
0621 internalId_.leadChEta_ = lLeadCh->eta();
0622 internalId_.leadChPhi_ = lLeadCh->phi();
0623 } else {
0624 internalId_.leadChPt_ = 0.0;
0625 internalId_.leadChEta_ = large_val;
0626 internalId_.leadChPhi_ = large_val;
0627 }
0628
0629 if (patjet != nullptr) {
0630 internalId_.nCharged_ = patjet->chargedMultiplicity();
0631 internalId_.nNeutrals_ = patjet->neutralMultiplicity();
0632 internalId_.chgEMfrac_ = patjet->chargedEmEnergy() / jet->energy();
0633 internalId_.neuEMfrac_ = patjet->neutralEmEnergy() / jet->energy();
0634 internalId_.chgHadrfrac_ = patjet->chargedHadronEnergy() / jet->energy();
0635 internalId_.neuHadrfrac_ = patjet->neutralHadronEnergy() / jet->energy();
0636 if (applyConstituentWeight)
0637 internalId_.nNeutrals_ = multNeut;
0638 } else {
0639 internalId_.nCharged_ = pfjet->chargedMultiplicity();
0640 internalId_.nNeutrals_ = pfjet->neutralMultiplicity();
0641 internalId_.chgEMfrac_ = pfjet->chargedEmEnergy() / jet->energy();
0642 internalId_.neuEMfrac_ = pfjet->neutralEmEnergy() / jet->energy();
0643 internalId_.chgHadrfrac_ = pfjet->chargedHadronEnergy() / jet->energy();
0644 internalId_.neuHadrfrac_ = pfjet->neutralHadronEnergy() / jet->energy();
0645 if (applyConstituentWeight)
0646 internalId_.nNeutrals_ = multNeut;
0647 }
0648 internalId_.nParticles_ = jet->nConstituents();
0649
0650
0651 float ddetaR_sum(0.0), ddphiR_sum(0.0), pull_tmp(0.0);
0652 for (unsigned k = 0; k < jet->numberOfSourceCandidatePtrs(); k++) {
0653 reco::CandidatePtr temp_pfJetConstituent = jet->sourceCandidatePtr(k);
0654
0655 const reco::Candidate* part = temp_pfJetConstituent.get();
0656
0657 float candWeight = 1.0;
0658
0659 if (applyConstituentWeight)
0660 candWeight = constituentWeights[jet->sourceCandidatePtr(k)];
0661
0662 float weight = candWeight * (part->pt()) * candWeight * (part->pt());
0663 float deta = part->eta() - jet->eta();
0664 float dphi = reco::deltaPhi(*part, *jet);
0665 float ddeta, ddphi, ddR;
0666 ddeta = deta - ave_deta;
0667 ddphi = dphi - ave_dphi;
0668 ddR = sqrt(ddeta * ddeta + ddphi * ddphi);
0669 ddetaR_sum += ddR * ddeta * weight;
0670 ddphiR_sum += ddR * ddphi * weight;
0671 }
0672 if (sumW2 > 0) {
0673 float ddetaR_ave = ddetaR_sum / sumW2;
0674 float ddphiR_ave = ddphiR_sum / sumW2;
0675 pull_tmp = sqrt(ddetaR_ave * ddetaR_ave + ddphiR_ave * ddphiR_ave);
0676 }
0677 internalId_.pull_ = pull_tmp;
0678
0679
0680 std::sort(frac.begin(), frac.end(), std::greater<float>());
0681 std::sort(fracCh.begin(), fracCh.end(), std::greater<float>());
0682 std::sort(fracEm.begin(), fracEm.end(), std::greater<float>());
0683 std::sort(fracNeut.begin(), fracNeut.end(), std::greater<float>());
0684 assign(frac, internalId_.leadFrac_, internalId_.secondFrac_, internalId_.thirdFrac_, internalId_.fourthFrac_);
0685 assign(
0686 fracCh, internalId_.leadChFrac_, internalId_.secondChFrac_, internalId_.thirdChFrac_, internalId_.fourthChFrac_);
0687 assign(
0688 fracEm, internalId_.leadEmFrac_, internalId_.secondEmFrac_, internalId_.thirdEmFrac_, internalId_.fourthEmFrac_);
0689 assign(fracNeut,
0690 internalId_.leadNeutFrac_,
0691 internalId_.secondNeutFrac_,
0692 internalId_.thirdNeutFrac_,
0693 internalId_.fourthNeutFrac_);
0694
0695 covMatrix(0, 0) /= sumPt2;
0696 covMatrix(0, 1) /= sumPt2;
0697 covMatrix(1, 1) /= sumPt2;
0698 covMatrix(1, 0) = covMatrix(0, 1);
0699 internalId_.etaW_ = sqrt(covMatrix(0, 0));
0700 internalId_.phiW_ = sqrt(covMatrix(1, 1));
0701 internalId_.jetW_ = 0.5 * (internalId_.etaW_ + internalId_.phiW_);
0702 TVectorD eigVals(2);
0703 eigVals = TMatrixDSymEigen(covMatrix).GetEigenValues();
0704 internalId_.majW_ = sqrt(std::abs(eigVals(0)));
0705 internalId_.minW_ = sqrt(std::abs(eigVals(1)));
0706 if (internalId_.majW_ < internalId_.minW_) {
0707 std::swap(internalId_.majW_, internalId_.minW_);
0708 }
0709
0710 internalId_.dRLeadCent_ = reco::deltaR(*jet, *lLead);
0711 if (lSecond != nullptr) {
0712 internalId_.dRLead2nd_ = reco::deltaR(*jet, *lSecond);
0713 }
0714 internalId_.dRMean_ /= jetPt;
0715 internalId_.dRMeanNeut_ /= jetPt;
0716 internalId_.dRMeanEm_ /= jetPt;
0717 internalId_.dRMeanCh_ /= jetPt;
0718 internalId_.dR2Mean_ /= sumPt2;
0719
0720 for (size_t ic = 0; ic < ncones; ++ic) {
0721 *coneFracs[ic] /= jetPt;
0722 *coneEmFracs[ic] /= jetPt;
0723 *coneNeutFracs[ic] /= jetPt;
0724 *coneChFracs[ic] /= jetPt;
0725 }
0726
0727 double ptMean = sumPt / internalId_.nParticles_;
0728 double ptRMS = 0;
0729 for (unsigned int i0 = 0; i0 < frac.size(); i0++) {
0730 ptRMS += (frac[i0] - ptMean) * (frac[i0] - ptMean);
0731 }
0732 ptRMS /= internalId_.nParticles_;
0733 ptRMS = sqrt(ptRMS);
0734
0735 internalId_.ptMean_ = ptMean;
0736 internalId_.ptRMS_ = ptRMS / jetPt;
0737 internalId_.pt2A_ = sqrt(internalId_.ptD_ / internalId_.nParticles_) / jetPt;
0738 internalId_.ptD_ = sqrt(internalId_.ptD_) / sumPt;
0739 internalId_.ptDCh_ = sqrt(internalId_.ptDCh_) / sumPtCh;
0740 internalId_.ptDNe_ = sqrt(internalId_.ptDNe_) / sumPtNe;
0741 internalId_.sumPt_ = sumPt;
0742 internalId_.sumChPt_ = sumPtCh;
0743 internalId_.sumNePt_ = sumPtNe;
0744
0745 internalId_.jetR_ = (lLead->pt()) * LeadcandWeight / sumPt;
0746 if (lLeadCh != nullptr) {
0747 internalId_.jetRchg_ = (lLeadCh->pt()) * LeadChcandWeight / sumPt;
0748 } else {
0749 internalId_.jetRchg_ = 0;
0750 }
0751
0752 internalId_.dRMatch_ = dRmin;
0753
0754 if (sumTkPt != 0.) {
0755 internalId_.beta_ /= sumTkPt;
0756 internalId_.betaStar_ /= sumTkPt;
0757 internalId_.betaClassic_ /= sumTkPt;
0758 internalId_.betaStarClassic_ /= sumTkPt;
0759 } else {
0760 assert(internalId_.beta_ == 0. && internalId_.betaStar_ == 0. && internalId_.betaClassic_ == 0. &&
0761 internalId_.betaStarClassic_ == 0.);
0762 }
0763
0764 if (cache_->runMvas()) {
0765 runMva();
0766 }
0767
0768 return PileupJetIdentifier(internalId_);
0769 }
0770
0771
0772 std::string PileupJetIdAlgo::dumpVariables() const {
0773 std::stringstream out;
0774 for (variables_list_t::const_iterator it = variables_.begin(); it != variables_.end(); ++it) {
0775 out << std::setw(15) << it->first << std::setw(3) << "=" << std::setw(5) << *it->second.first << " ("
0776 << std::setw(5) << it->second.second << ")" << std::endl;
0777 }
0778 return out.str();
0779 }
0780
0781
0782 void PileupJetIdAlgo::resetVariables() {
0783 internalId_.idFlag_ = 0;
0784 for (variables_list_t::iterator it = variables_.begin(); it != variables_.end(); ++it) {
0785 *it->second.first = it->second.second;
0786 }
0787 }
0788
0789
0790 #define INIT_VARIABLE(NAME, TMVANAME, VAL) \
0791 internalId_.NAME##_ = VAL; \
0792 variables_[#NAME] = std::make_pair(&internalId_.NAME##_, VAL);
0793
0794
0795 void PileupJetIdAlgo::initVariables() {
0796 internalId_.idFlag_ = 0;
0797 INIT_VARIABLE(mva, "", -100.);
0798
0799
0800 INIT_VARIABLE(jetPt, "", 0.);
0801 INIT_VARIABLE(jetEta, "", large_val);
0802 INIT_VARIABLE(jetPhi, "jsphi_1", large_val);
0803 INIT_VARIABLE(jetM, "jm_1", 0.);
0804
0805 INIT_VARIABLE(nCharged, "", 0.);
0806 INIT_VARIABLE(nNeutrals, "", 0.);
0807
0808 INIT_VARIABLE(chgEMfrac, "", 0.);
0809 INIT_VARIABLE(neuEMfrac, "", 0.);
0810 INIT_VARIABLE(chgHadrfrac, "", 0.);
0811 INIT_VARIABLE(neuHadrfrac, "", 0.);
0812
0813 INIT_VARIABLE(d0, "jd0_1", -1000.);
0814 INIT_VARIABLE(dZ, "jdZ_1", -1000.);
0815
0816 INIT_VARIABLE(nParticles, "", 0.);
0817
0818 INIT_VARIABLE(leadPt, "lpt_1", 0.);
0819 INIT_VARIABLE(leadEta, "leta_1", large_val);
0820 INIT_VARIABLE(leadPhi, "lphi_1", large_val);
0821 INIT_VARIABLE(secondPt, "spt_1", 0.);
0822 INIT_VARIABLE(secondEta, "seta_1", large_val);
0823 INIT_VARIABLE(secondPhi, "sphi_1", large_val);
0824 INIT_VARIABLE(leadNeutPt, "lnept_1", 0.);
0825 INIT_VARIABLE(leadNeutEta, "lneeta_1", large_val);
0826 INIT_VARIABLE(leadNeutPhi, "lnephi_1", large_val);
0827 INIT_VARIABLE(leadEmPt, "lempt_1", 0.);
0828 INIT_VARIABLE(leadEmEta, "lemeta_1", large_val);
0829 INIT_VARIABLE(leadEmPhi, "lemphi_1", large_val);
0830 INIT_VARIABLE(leadChPt, "lchpt_1", 0.);
0831 INIT_VARIABLE(leadChEta, "lcheta_1", large_val);
0832 INIT_VARIABLE(leadChPhi, "lchphi_1", large_val);
0833 INIT_VARIABLE(leadFrac, "lLfr_1", 0.);
0834
0835 INIT_VARIABLE(dRLeadCent, "drlc_1", 0.);
0836 INIT_VARIABLE(dRLead2nd, "drls_1", 0.);
0837 INIT_VARIABLE(dRMean, "drm_1", 0.);
0838 INIT_VARIABLE(dRMean, "", 0.);
0839 INIT_VARIABLE(pull, "", 0.);
0840 INIT_VARIABLE(dRMeanNeut, "drmne_1", 0.);
0841 INIT_VARIABLE(dRMeanEm, "drem_1", 0.);
0842 INIT_VARIABLE(dRMeanCh, "drch_1", 0.);
0843 INIT_VARIABLE(dR2Mean, "", 0.);
0844
0845 INIT_VARIABLE(ptD, "", 0.);
0846 INIT_VARIABLE(ptMean, "", 0.);
0847 INIT_VARIABLE(ptRMS, "", 0.);
0848 INIT_VARIABLE(pt2A, "", 0.);
0849 INIT_VARIABLE(ptDCh, "", 0.);
0850 INIT_VARIABLE(ptDNe, "", 0.);
0851 INIT_VARIABLE(sumPt, "", 0.);
0852 INIT_VARIABLE(sumChPt, "", 0.);
0853 INIT_VARIABLE(sumNePt, "", 0.);
0854
0855 INIT_VARIABLE(secondFrac, "", 0.);
0856 INIT_VARIABLE(thirdFrac, "", 0.);
0857 INIT_VARIABLE(fourthFrac, "", 0.);
0858
0859 INIT_VARIABLE(leadChFrac, "", 0.);
0860 INIT_VARIABLE(secondChFrac, "", 0.);
0861 INIT_VARIABLE(thirdChFrac, "", 0.);
0862 INIT_VARIABLE(fourthChFrac, "", 0.);
0863
0864 INIT_VARIABLE(leadNeutFrac, "", 0.);
0865 INIT_VARIABLE(secondNeutFrac, "", 0.);
0866 INIT_VARIABLE(thirdNeutFrac, "", 0.);
0867 INIT_VARIABLE(fourthNeutFrac, "", 0.);
0868
0869 INIT_VARIABLE(leadEmFrac, "", 0.);
0870 INIT_VARIABLE(secondEmFrac, "", 0.);
0871 INIT_VARIABLE(thirdEmFrac, "", 0.);
0872 INIT_VARIABLE(fourthEmFrac, "", 0.);
0873
0874 INIT_VARIABLE(jetW, "", 1.);
0875 INIT_VARIABLE(etaW, "", 1.);
0876 INIT_VARIABLE(phiW, "", 1.);
0877
0878 INIT_VARIABLE(majW, "", 1.);
0879 INIT_VARIABLE(minW, "", 1.);
0880
0881 INIT_VARIABLE(frac01, "", 0.);
0882 INIT_VARIABLE(frac02, "", 0.);
0883 INIT_VARIABLE(frac03, "", 0.);
0884 INIT_VARIABLE(frac04, "", 0.);
0885 INIT_VARIABLE(frac05, "", 0.);
0886 INIT_VARIABLE(frac06, "", 0.);
0887 INIT_VARIABLE(frac07, "", 0.);
0888
0889 INIT_VARIABLE(chFrac01, "", 0.);
0890 INIT_VARIABLE(chFrac02, "", 0.);
0891 INIT_VARIABLE(chFrac03, "", 0.);
0892 INIT_VARIABLE(chFrac04, "", 0.);
0893 INIT_VARIABLE(chFrac05, "", 0.);
0894 INIT_VARIABLE(chFrac06, "", 0.);
0895 INIT_VARIABLE(chFrac07, "", 0.);
0896
0897 INIT_VARIABLE(neutFrac01, "", 0.);
0898 INIT_VARIABLE(neutFrac02, "", 0.);
0899 INIT_VARIABLE(neutFrac03, "", 0.);
0900 INIT_VARIABLE(neutFrac04, "", 0.);
0901 INIT_VARIABLE(neutFrac05, "", 0.);
0902 INIT_VARIABLE(neutFrac06, "", 0.);
0903 INIT_VARIABLE(neutFrac07, "", 0.);
0904
0905 INIT_VARIABLE(emFrac01, "", 0.);
0906 INIT_VARIABLE(emFrac02, "", 0.);
0907 INIT_VARIABLE(emFrac03, "", 0.);
0908 INIT_VARIABLE(emFrac04, "", 0.);
0909 INIT_VARIABLE(emFrac05, "", 0.);
0910 INIT_VARIABLE(emFrac06, "", 0.);
0911 INIT_VARIABLE(emFrac07, "", 0.);
0912
0913 INIT_VARIABLE(beta, "", 0.);
0914 INIT_VARIABLE(betaStar, "", 0.);
0915 INIT_VARIABLE(betaClassic, "", 0.);
0916 INIT_VARIABLE(betaStarClassic, "", 0.);
0917
0918 INIT_VARIABLE(nvtx, "", 0.);
0919 INIT_VARIABLE(rho, "", 0.);
0920 INIT_VARIABLE(nTrueInt, "", 0.);
0921
0922 INIT_VARIABLE(jetR, "", 0.);
0923 INIT_VARIABLE(jetRchg, "", 0.);
0924 INIT_VARIABLE(dRMatch, "", 0.);
0925 }
0926
0927 #undef INIT_VARIABLE
