Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 09:59:53

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   // initialize all variables to 0
0269   resetVariables();
0270 
0271   // loop over constituents, accumulate sums and find leading candidates
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.) {  // if this is a pat jet and no jec has been passed take the jec from the object
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;  // use uncorrected pt for shape variables
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_);  // use corrected pt for jet kinematics
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     const pat::PackedCandidate* lPack = dynamic_cast<const pat::PackedCandidate*>(icand);
0342     const reco::PFCandidate* lPF = dynamic_cast<const reco::PFCandidate*>(icand);
0343     bool isPacked = true;
0344     if (lPack == nullptr) {
0345       isPacked = false;
0346     }
0347 
0348     float candWeight = 1.0;
0349     if (applyConstituentWeight) {  // PUPPI Jet weight should be pulled up from valuemap, not packed candidate
0350       candWeight = constituentWeights[jet->sourceCandidatePtr(i)];
0351     }
0352     float candPt = (icand->pt()) * candWeight;
0353     float candPtFrac = candPt / jetPt;
0354     float candDr = reco::deltaR(*icand, *jet);
0355     float candDeta = icand->eta() - jet->eta();
0356     float candDphi = reco::deltaPhi(*icand, *jet);
0357     float candPtDr = candPt * candDr;
0358     size_t icone = std::lower_bound(&cones[0], &cones[ncones], candDr) - &cones[0];
0359 
0360     if (candDr < dRmin)
0361       dRmin = candDr;
0362 
0363     // // all particles; PUPPI weights multiplied to leading and subleading constituent if it is for PUPPI
0364     if (lLead == nullptr || candPt > (lLead->pt()) * LeadcandWeight) {
0365       lSecond = lLead;
0366       SecondcandWeight = LeadcandWeight;
0367       lLead = icand;
0368       if (applyConstituentWeight) {
0369         LeadcandWeight = constituentWeights[jet->sourceCandidatePtr(i)];
0370       }
0371     } else if ((lSecond == nullptr || candPt > (lSecond->pt()) * SecondcandWeight) &&
0372                (candPt < (lLead->pt()) * LeadcandWeight)) {
0373       lSecond = icand;
0374       if (applyConstituentWeight) {
0375         SecondcandWeight = constituentWeights[jet->sourceCandidatePtr(i)];
0376       }
0377     }
0378 
0379     // // average shapes
0380     internalId_.dRMean_ += candPtDr;
0381     internalId_.dR2Mean_ += candPtDr * candPtDr;
0382     covMatrix(0, 0) += candPt * candPt * candDeta * candDeta;
0383     covMatrix(0, 1) += candPt * candPt * candDeta * candDphi;
0384     covMatrix(1, 1) += candPt * candPt * candDphi * candDphi;
0385     internalId_.ptD_ += candPt * candPt;
0386     sumPt += candPt;
0387     sumPt2 += candPt * candPt;
0388 
0389     // single most energetic candiates and jet shape profiles
0390     frac.push_back(candPtFrac);
0391 
0392     if (icone < ncones) {
0393       *coneFracs[icone] += candPt;
0394     }
0395 
0396     // neutrals Neutral hadrons
0397     if (abs(icand->pdgId()) == 130) {
0398       if (lLeadNeut == nullptr || candPt > (lLeadNeut->pt()) * LeadNeutcandWeight) {
0399         lLeadNeut = icand;
0400         if (applyConstituentWeight) {
0401           LeadNeutcandWeight = constituentWeights[jet->sourceCandidatePtr(i)];
0402         }
0403       }
0404 
0405       internalId_.dRMeanNeut_ += candPtDr;
0406       fracNeut.push_back(candPtFrac);
0407       if (icone < ncones) {
0408         *coneNeutFracs[icone] += candPt;
0409       }
0410       internalId_.ptDNe_ += candPt * candPt;
0411       sumPtNe += candPt;
0412       multNeut += candWeight;
0413     }
0414 
0415     // EM candidated photon
0416     if (icand->pdgId() == 22) {
0417       if (lLeadEm == nullptr || candPt > (lLeadEm->pt()) * LeadEmcandWeight) {
0418         lLeadEm = icand;
0419         if (applyConstituentWeight) {
0420           LeadEmcandWeight = constituentWeights[jet->sourceCandidatePtr(i)];
0421         }
0422       }
0423       internalId_.dRMeanEm_ += candPtDr;
0424       fracEm.push_back(candPtFrac);
0425       if (icone < ncones) {
0426         *coneEmFracs[icone] += candPt;
0427       }
0428       internalId_.ptDNe_ += candPt * candPt;
0429       sumPtNe += candPt;
0430       multNeut += candWeight;
0431     }
0432     // hadrons and EM in HF
0433     if ((abs(icand->pdgId()) == 1) || (abs(icand->pdgId()) == 2))
0434       multNeut += candWeight;
0435 
0436     // Charged  particles
0437     if (icand->charge() != 0) {
0438       if (lLeadCh == nullptr || candPt > (lLeadCh->pt()) * LeadChcandWeight) {
0439         lLeadCh = icand;
0440         if (applyConstituentWeight) {
0441           LeadChcandWeight = constituentWeights[jet->sourceCandidatePtr(i)];
0442         }
0443         const reco::Track* pfTrk = icand->bestTrack();
0444         if (lPF && std::abs(icand->pdgId()) == 13 && pfTrk == nullptr) {
0445           reco::MuonRef lmuRef = lPF->muonRef();
0446           if (lmuRef.isNonnull()) {
0447             const reco::Muon& lmu = *lmuRef.get();
0448             pfTrk = lmu.bestTrack();
0449             edm::LogWarning("BadMuon")
0450                 << "Found a PFCandidate muon without a trackRef: falling back to Muon::bestTrack ";
0451           }
0452         }
0453         if (pfTrk == nullptr) {  //protection against empty pointers for the miniAOD case
0454           //To handle the electron case
0455           if (isPacked) {
0456             internalId_.d0_ = std::abs(lPack->dxy(vtx->position()));
0457             internalId_.dZ_ = std::abs(lPack->dz(vtx->position()));
0458           } else if (lPF != nullptr) {
0459             pfTrk = (lPF->trackRef().get() == nullptr) ? lPF->gsfTrackRef().get() : lPF->trackRef().get();
0460             internalId_.d0_ = std::abs(pfTrk->dxy(vtx->position()));
0461             internalId_.dZ_ = std::abs(pfTrk->dz(vtx->position()));
0462           }
0463         } else {
0464           internalId_.d0_ = std::abs(pfTrk->dxy(vtx->position()));
0465           internalId_.dZ_ = std::abs(pfTrk->dz(vtx->position()));
0466         }
0467       }
0468       internalId_.dRMeanCh_ += candPtDr;
0469       internalId_.ptDCh_ += candPt * candPt;
0470       fracCh.push_back(candPtFrac);
0471       if (icone < ncones) {
0472         *coneChFracs[icone] += candPt;
0473       }
0474       sumPtCh += candPt;
0475     }
0476     // // beta and betastar
0477     if (icand->charge() != 0) {
0478       if (!isPacked) {
0479         if (lPF->trackRef().isNonnull()) {
0480           float tkpt = candPt;
0481           sumTkPt += tkpt;
0482           // 'classic' beta definition based on track-vertex association
0483           bool inVtx0 = vtx->trackWeight(lPF->trackRef()) > 0;
0484 
0485           bool inAnyOther = false;
0486           // alternative beta definition based on track-vertex distance of closest approach
0487           double dZ0 = std::abs(lPF->trackRef()->dz(vtx->position()));
0488           double dZ = dZ0;
0489           for (reco::VertexCollection::const_iterator vi = allvtx.begin(); vi != allvtx.end(); ++vi) {
0490             const reco::Vertex& iv = *vi;
0491             if (iv.isFake() || iv.ndof() < 4) {
0492               continue;
0493             }
0494             // the primary vertex may have been copied by the user: check identity by position
0495             bool isVtx0 = (iv.position() - vtx->position()).r() < 0.02;
0496             // 'classic' beta definition: check if the track is associated with
0497             // any vertex other than the primary one
0498             if (!isVtx0 && !inAnyOther) {
0499               inAnyOther = vtx->trackWeight(lPF->trackRef()) <= 0;
0500             }
0501             // alternative beta: find closest vertex to the track
0502             dZ = std::min(dZ, std::abs(lPF->trackRef()->dz(iv.position())));
0503           }
0504           // classic beta/betaStar
0505           if (inVtx0 && !inAnyOther) {
0506             internalId_.betaClassic_ += tkpt;
0507           } else if (!inVtx0 && inAnyOther) {
0508             internalId_.betaStarClassic_ += tkpt;
0509           }
0510           // alternative beta/betaStar
0511           if (dZ0 < 0.2) {
0512             internalId_.beta_ += tkpt;
0513           } else if (dZ < 0.2) {
0514             internalId_.betaStar_ += tkpt;
0515           }
0516         }
0517       } else {
0518         float tkpt = candPt;
0519         sumTkPt += tkpt;
0520         bool inVtx0 = false;
0521         bool inVtxOther = false;
0522         double dZ0 = 9999.;
0523         double dZ_tmp = 9999.;
0524         for (unsigned vtx_i = 0; vtx_i < allvtx.size(); vtx_i++) {
0525           const auto& iv = allvtx[vtx_i];
0526 
0527           if (iv.isFake())
0528             continue;
0529 
0530           // Match to vertex in case of copy as above
0531           bool isVtx0 = (iv.position() - vtx->position()).r() < 0.02;
0532 
0533           if (isVtx0) {
0534             if (lPack->fromPV(vtx_i) == pat::PackedCandidate::PVUsedInFit)
0535               inVtx0 = true;
0536             if (lPack->fromPV(vtx_i) == 0)
0537               inVtxOther = true;
0538             dZ0 = lPack->dz(iv.position());
0539           }
0540 
0541           if (fabs(lPack->dz(iv.position())) < fabs(dZ_tmp)) {
0542             dZ_tmp = lPack->dz(iv.position());
0543           }
0544         }
0545         if (inVtx0) {
0546           internalId_.betaClassic_ += tkpt;
0547         } else if (inVtxOther) {
0548           internalId_.betaStarClassic_ += tkpt;
0549         }
0550         if (fabs(dZ0) < 0.2) {
0551           internalId_.beta_ += tkpt;
0552         } else if (fabs(dZ_tmp) < 0.2) {
0553           internalId_.betaStar_ += tkpt;
0554         }
0555       }
0556     }
0557 
0558     // trailing candidate
0559     if (lTrail == nullptr || candPt < (lTrail->pt()) * TrailcandWeight) {
0560       lTrail = icand;
0561       if (applyConstituentWeight) {
0562         TrailcandWeight = constituentWeights[jet->sourceCandidatePtr(i)];
0563       }
0564     }
0565 
0566     // average for pull variable
0567 
0568     float weight2 = candPt * candPt;
0569     sumW2 += weight2;
0570     float deta = icand->eta() - jet->eta();
0571     float dphi = reco::deltaPhi(*icand, *jet);
0572     sum_deta += deta * weight2;
0573     sum_dphi += dphi * weight2;
0574   }
0575   if (sumW2 > 0) {
0576     ave_deta = sum_deta / sumW2;
0577     ave_dphi = sum_dphi / sumW2;
0578   }
0579 
0580   // // Finalize all variables
0581   // Most of Below values are not used for puID variable generation at the moment, except lLeadCh Pt for JetRchg, so I assign that zero if there is no charged constituent.
0582 
0583   assert(!(lLead == nullptr));
0584   internalId_.leadPt_ = lLead->pt() * LeadcandWeight;
0585   internalId_.leadEta_ = lLead->eta();
0586   internalId_.leadPhi_ = lLead->phi();
0587 
0588   if (lSecond != nullptr) {
0589     internalId_.secondPt_ = lSecond->pt() * SecondcandWeight;
0590     internalId_.secondEta_ = lSecond->eta();
0591     internalId_.secondPhi_ = lSecond->phi();
0592   } else {
0593     internalId_.secondPt_ = 0.0;
0594     internalId_.secondEta_ = large_val;
0595     internalId_.secondPhi_ = large_val;
0596   }
0597 
0598   if (lLeadNeut != nullptr) {
0599     internalId_.leadNeutPt_ = lLeadNeut->pt() * LeadNeutcandWeight;
0600     internalId_.leadNeutEta_ = lLeadNeut->eta();
0601     internalId_.leadNeutPhi_ = lLeadNeut->phi();
0602   } else {
0603     internalId_.leadNeutPt_ = 0.0;
0604     internalId_.leadNeutEta_ = large_val;
0605     internalId_.leadNeutPhi_ = large_val;
0606   }
0607 
0608   if (lLeadEm != nullptr) {
0609     internalId_.leadEmPt_ = lLeadEm->pt() * LeadEmcandWeight;
0610     internalId_.leadEmEta_ = lLeadEm->eta();
0611     internalId_.leadEmPhi_ = lLeadEm->phi();
0612   } else {
0613     internalId_.leadEmPt_ = 0.0;
0614     internalId_.leadEmEta_ = large_val;
0615     internalId_.leadEmPhi_ = large_val;
0616   }
0617 
0618   if (lLeadCh != nullptr) {
0619     internalId_.leadChPt_ = lLeadCh->pt() * LeadChcandWeight;
0620     internalId_.leadChEta_ = lLeadCh->eta();
0621     internalId_.leadChPhi_ = lLeadCh->phi();
0622   } else {
0623     internalId_.leadChPt_ = 0.0;
0624     internalId_.leadChEta_ = large_val;
0625     internalId_.leadChPhi_ = large_val;
0626   }
0627 
0628   if (patjet != nullptr) {  // to enable running on MiniAOD slimmedJets
0629     internalId_.nCharged_ = patjet->chargedMultiplicity();
0630     internalId_.nNeutrals_ = patjet->neutralMultiplicity();
0631     internalId_.chgEMfrac_ = patjet->chargedEmEnergy() / jet->energy();
0632     internalId_.neuEMfrac_ = patjet->neutralEmEnergy() / jet->energy();
0633     internalId_.chgHadrfrac_ = patjet->chargedHadronEnergy() / jet->energy();
0634     internalId_.neuHadrfrac_ = patjet->neutralHadronEnergy() / jet->energy();
0635     if (applyConstituentWeight)
0636       internalId_.nNeutrals_ = multNeut;
0637   } else {
0638     internalId_.nCharged_ = pfjet->chargedMultiplicity();
0639     internalId_.nNeutrals_ = pfjet->neutralMultiplicity();
0640     internalId_.chgEMfrac_ = pfjet->chargedEmEnergy() / jet->energy();
0641     internalId_.neuEMfrac_ = pfjet->neutralEmEnergy() / jet->energy();
0642     internalId_.chgHadrfrac_ = pfjet->chargedHadronEnergy() / jet->energy();
0643     internalId_.neuHadrfrac_ = pfjet->neutralHadronEnergy() / jet->energy();
0644     if (applyConstituentWeight)
0645       internalId_.nNeutrals_ = multNeut;
0646   }
0647   internalId_.nParticles_ = jet->nConstituents();
0648 
0649   ///////////////////////pull variable///////////////////////////////////
0650   float ddetaR_sum(0.0), ddphiR_sum(0.0), pull_tmp(0.0);
0651   for (unsigned k = 0; k < jet->numberOfSourceCandidatePtrs(); k++) {
0652     reco::CandidatePtr temp_pfJetConstituent = jet->sourceCandidatePtr(k);
0653     //    reco::CandidatePtr temp_weightpfJetConstituent = jet->sourceCandidatePtr(k);
0654     const reco::Candidate* part = temp_pfJetConstituent.get();
0655 
0656     float candWeight = 1.0;
0657 
0658     if (applyConstituentWeight)
0659       candWeight = constituentWeights[jet->sourceCandidatePtr(k)];
0660 
0661     float weight = candWeight * (part->pt()) * candWeight * (part->pt());
0662     float deta = part->eta() - jet->eta();
0663     float dphi = reco::deltaPhi(*part, *jet);
0664     float ddeta, ddphi, ddR;
0665     ddeta = deta - ave_deta;
0666     ddphi = dphi - ave_dphi;
0667     ddR = sqrt(ddeta * ddeta + ddphi * ddphi);
0668     ddetaR_sum += ddR * ddeta * weight;
0669     ddphiR_sum += ddR * ddphi * weight;
0670   }
0671   if (sumW2 > 0) {
0672     float ddetaR_ave = ddetaR_sum / sumW2;
0673     float ddphiR_ave = ddphiR_sum / sumW2;
0674     pull_tmp = sqrt(ddetaR_ave * ddetaR_ave + ddphiR_ave * ddphiR_ave);
0675   }
0676   internalId_.pull_ = pull_tmp;
0677   ///////////////////////////////////////////////////////////////////////
0678 
0679   std::sort(frac.begin(), frac.end(), std::greater<float>());
0680   std::sort(fracCh.begin(), fracCh.end(), std::greater<float>());
0681   std::sort(fracEm.begin(), fracEm.end(), std::greater<float>());
0682   std::sort(fracNeut.begin(), fracNeut.end(), std::greater<float>());
0683   assign(frac, internalId_.leadFrac_, internalId_.secondFrac_, internalId_.thirdFrac_, internalId_.fourthFrac_);
0684   assign(
0685       fracCh, internalId_.leadChFrac_, internalId_.secondChFrac_, internalId_.thirdChFrac_, internalId_.fourthChFrac_);
0686   assign(
0687       fracEm, internalId_.leadEmFrac_, internalId_.secondEmFrac_, internalId_.thirdEmFrac_, internalId_.fourthEmFrac_);
0688   assign(fracNeut,
0689          internalId_.leadNeutFrac_,
0690          internalId_.secondNeutFrac_,
0691          internalId_.thirdNeutFrac_,
0692          internalId_.fourthNeutFrac_);
0693 
0694   covMatrix(0, 0) /= sumPt2;
0695   covMatrix(0, 1) /= sumPt2;
0696   covMatrix(1, 1) /= sumPt2;
0697   covMatrix(1, 0) = covMatrix(0, 1);
0698   internalId_.etaW_ = sqrt(covMatrix(0, 0));
0699   internalId_.phiW_ = sqrt(covMatrix(1, 1));
0700   internalId_.jetW_ = 0.5 * (internalId_.etaW_ + internalId_.phiW_);
0701   TVectorD eigVals(2);
0702   eigVals = TMatrixDSymEigen(covMatrix).GetEigenValues();
0703   internalId_.majW_ = sqrt(std::abs(eigVals(0)));
0704   internalId_.minW_ = sqrt(std::abs(eigVals(1)));
0705   if (internalId_.majW_ < internalId_.minW_) {
0706     std::swap(internalId_.majW_, internalId_.minW_);
0707   }
0708 
0709   internalId_.dRLeadCent_ = reco::deltaR(*jet, *lLead);
0710   if (lSecond != nullptr) {
0711     internalId_.dRLead2nd_ = reco::deltaR(*jet, *lSecond);
0712   }
0713   internalId_.dRMean_ /= jetPt;
0714   internalId_.dRMeanNeut_ /= jetPt;
0715   internalId_.dRMeanEm_ /= jetPt;
0716   internalId_.dRMeanCh_ /= jetPt;
0717   internalId_.dR2Mean_ /= sumPt2;
0718 
0719   for (size_t ic = 0; ic < ncones; ++ic) {
0720     *coneFracs[ic] /= jetPt;
0721     *coneEmFracs[ic] /= jetPt;
0722     *coneNeutFracs[ic] /= jetPt;
0723     *coneChFracs[ic] /= jetPt;
0724   }
0725   //http://jets.physics.harvard.edu/qvg/
0726   double ptMean = sumPt / internalId_.nParticles_;
0727   double ptRMS = 0;
0728   for (unsigned int i0 = 0; i0 < frac.size(); i0++) {
0729     ptRMS += (frac[i0] - ptMean) * (frac[i0] - ptMean);
0730   }
0731   ptRMS /= internalId_.nParticles_;
0732   ptRMS = sqrt(ptRMS);
0733 
0734   internalId_.ptMean_ = ptMean;
0735   internalId_.ptRMS_ = ptRMS / jetPt;
0736   internalId_.pt2A_ = sqrt(internalId_.ptD_ / internalId_.nParticles_) / jetPt;
0737   internalId_.ptD_ = sqrt(internalId_.ptD_) / sumPt;
0738   internalId_.ptDCh_ = sqrt(internalId_.ptDCh_) / sumPtCh;
0739   internalId_.ptDNe_ = sqrt(internalId_.ptDNe_) / sumPtNe;
0740   internalId_.sumPt_ = sumPt;
0741   internalId_.sumChPt_ = sumPtCh;
0742   internalId_.sumNePt_ = sumPtNe;
0743 
0744   internalId_.jetR_ = (lLead->pt()) * LeadcandWeight / sumPt;
0745   if (lLeadCh != nullptr) {
0746     internalId_.jetRchg_ = (lLeadCh->pt()) * LeadChcandWeight / sumPt;
0747   } else {
0748     internalId_.jetRchg_ = 0;
0749   }
0750 
0751   internalId_.dRMatch_ = dRmin;
0752 
0753   if (sumTkPt != 0.) {
0754     internalId_.beta_ /= sumTkPt;
0755     internalId_.betaStar_ /= sumTkPt;
0756     internalId_.betaClassic_ /= sumTkPt;
0757     internalId_.betaStarClassic_ /= sumTkPt;
0758   } else {
0759     assert(internalId_.beta_ == 0. && internalId_.betaStar_ == 0. && internalId_.betaClassic_ == 0. &&
0760            internalId_.betaStarClassic_ == 0.);
0761   }
0762 
0763   if (cache_->runMvas()) {
0764     runMva();
0765   }
0766 
0767   return PileupJetIdentifier(internalId_);
0768 }
0769 
0770 // ------------------------------------------------------------------------------------------
0771 std::string PileupJetIdAlgo::dumpVariables() const {
0772   std::stringstream out;
0773   for (variables_list_t::const_iterator it = variables_.begin(); it != variables_.end(); ++it) {
0774     out << std::setw(15) << it->first << std::setw(3) << "=" << std::setw(5) << *it->second.first << " ("
0775         << std::setw(5) << it->second.second << ")" << std::endl;
0776   }
0777   return out.str();
0778 }
0779 
0780 // ------------------------------------------------------------------------------------------
0781 void PileupJetIdAlgo::resetVariables() {
0782   internalId_.idFlag_ = 0;
0783   for (variables_list_t::iterator it = variables_.begin(); it != variables_.end(); ++it) {
0784     *it->second.first = it->second.second;
0785   }
0786 }
0787 
0788 // ------------------------------------------------------------------------------------------
0789 #define INIT_VARIABLE(NAME, TMVANAME, VAL) \
0790   internalId_.NAME##_ = VAL;               \
0791   variables_[#NAME] = std::make_pair(&internalId_.NAME##_, VAL);
0792 
0793 // ------------------------------------------------------------------------------------------
0794 void PileupJetIdAlgo::initVariables() {
0795   internalId_.idFlag_ = 0;
0796   INIT_VARIABLE(mva, "", -100.);
0797   //INIT_VARIABLE(jetPt      , "jspt_1", 0.);
0798   //INIT_VARIABLE(jetEta     , "jseta_1", large_val);
0799   INIT_VARIABLE(jetPt, "", 0.);
0800   INIT_VARIABLE(jetEta, "", large_val);
0801   INIT_VARIABLE(jetPhi, "jsphi_1", large_val);
0802   INIT_VARIABLE(jetM, "jm_1", 0.);
0803 
0804   INIT_VARIABLE(nCharged, "", 0.);
0805   INIT_VARIABLE(nNeutrals, "", 0.);
0806 
0807   INIT_VARIABLE(chgEMfrac, "", 0.);
0808   INIT_VARIABLE(neuEMfrac, "", 0.);
0809   INIT_VARIABLE(chgHadrfrac, "", 0.);
0810   INIT_VARIABLE(neuHadrfrac, "", 0.);
0811 
0812   INIT_VARIABLE(d0, "jd0_1", -1000.);
0813   INIT_VARIABLE(dZ, "jdZ_1", -1000.);
0814   //INIT_VARIABLE(nParticles , "npart_1"  , 0.);
0815   INIT_VARIABLE(nParticles, "", 0.);
0816 
0817   INIT_VARIABLE(leadPt, "lpt_1", 0.);
0818   INIT_VARIABLE(leadEta, "leta_1", large_val);
0819   INIT_VARIABLE(leadPhi, "lphi_1", large_val);
0820   INIT_VARIABLE(secondPt, "spt_1", 0.);
0821   INIT_VARIABLE(secondEta, "seta_1", large_val);
0822   INIT_VARIABLE(secondPhi, "sphi_1", large_val);
0823   INIT_VARIABLE(leadNeutPt, "lnept_1", 0.);
0824   INIT_VARIABLE(leadNeutEta, "lneeta_1", large_val);
0825   INIT_VARIABLE(leadNeutPhi, "lnephi_1", large_val);
0826   INIT_VARIABLE(leadEmPt, "lempt_1", 0.);
0827   INIT_VARIABLE(leadEmEta, "lemeta_1", large_val);
0828   INIT_VARIABLE(leadEmPhi, "lemphi_1", large_val);
0829   INIT_VARIABLE(leadChPt, "lchpt_1", 0.);
0830   INIT_VARIABLE(leadChEta, "lcheta_1", large_val);
0831   INIT_VARIABLE(leadChPhi, "lchphi_1", large_val);
0832   INIT_VARIABLE(leadFrac, "lLfr_1", 0.);
0833 
0834   INIT_VARIABLE(dRLeadCent, "drlc_1", 0.);
0835   INIT_VARIABLE(dRLead2nd, "drls_1", 0.);
0836   INIT_VARIABLE(dRMean, "drm_1", 0.);
0837   INIT_VARIABLE(dRMean, "", 0.);
0838   INIT_VARIABLE(pull, "", 0.);
0839   INIT_VARIABLE(dRMeanNeut, "drmne_1", 0.);
0840   INIT_VARIABLE(dRMeanEm, "drem_1", 0.);
0841   INIT_VARIABLE(dRMeanCh, "drch_1", 0.);
0842   INIT_VARIABLE(dR2Mean, "", 0.);
0843 
0844   INIT_VARIABLE(ptD, "", 0.);
0845   INIT_VARIABLE(ptMean, "", 0.);
0846   INIT_VARIABLE(ptRMS, "", 0.);
0847   INIT_VARIABLE(pt2A, "", 0.);
0848   INIT_VARIABLE(ptDCh, "", 0.);
0849   INIT_VARIABLE(ptDNe, "", 0.);
0850   INIT_VARIABLE(sumPt, "", 0.);
0851   INIT_VARIABLE(sumChPt, "", 0.);
0852   INIT_VARIABLE(sumNePt, "", 0.);
0853 
0854   INIT_VARIABLE(secondFrac, "", 0.);
0855   INIT_VARIABLE(thirdFrac, "", 0.);
0856   INIT_VARIABLE(fourthFrac, "", 0.);
0857 
0858   INIT_VARIABLE(leadChFrac, "", 0.);
0859   INIT_VARIABLE(secondChFrac, "", 0.);
0860   INIT_VARIABLE(thirdChFrac, "", 0.);
0861   INIT_VARIABLE(fourthChFrac, "", 0.);
0862 
0863   INIT_VARIABLE(leadNeutFrac, "", 0.);
0864   INIT_VARIABLE(secondNeutFrac, "", 0.);
0865   INIT_VARIABLE(thirdNeutFrac, "", 0.);
0866   INIT_VARIABLE(fourthNeutFrac, "", 0.);
0867 
0868   INIT_VARIABLE(leadEmFrac, "", 0.);
0869   INIT_VARIABLE(secondEmFrac, "", 0.);
0870   INIT_VARIABLE(thirdEmFrac, "", 0.);
0871   INIT_VARIABLE(fourthEmFrac, "", 0.);
0872 
0873   INIT_VARIABLE(jetW, "", 1.);
0874   INIT_VARIABLE(etaW, "", 1.);
0875   INIT_VARIABLE(phiW, "", 1.);
0876 
0877   INIT_VARIABLE(majW, "", 1.);
0878   INIT_VARIABLE(minW, "", 1.);
0879 
0880   INIT_VARIABLE(frac01, "", 0.);
0881   INIT_VARIABLE(frac02, "", 0.);
0882   INIT_VARIABLE(frac03, "", 0.);
0883   INIT_VARIABLE(frac04, "", 0.);
0884   INIT_VARIABLE(frac05, "", 0.);
0885   INIT_VARIABLE(frac06, "", 0.);
0886   INIT_VARIABLE(frac07, "", 0.);
0887 
0888   INIT_VARIABLE(chFrac01, "", 0.);
0889   INIT_VARIABLE(chFrac02, "", 0.);
0890   INIT_VARIABLE(chFrac03, "", 0.);
0891   INIT_VARIABLE(chFrac04, "", 0.);
0892   INIT_VARIABLE(chFrac05, "", 0.);
0893   INIT_VARIABLE(chFrac06, "", 0.);
0894   INIT_VARIABLE(chFrac07, "", 0.);
0895 
0896   INIT_VARIABLE(neutFrac01, "", 0.);
0897   INIT_VARIABLE(neutFrac02, "", 0.);
0898   INIT_VARIABLE(neutFrac03, "", 0.);
0899   INIT_VARIABLE(neutFrac04, "", 0.);
0900   INIT_VARIABLE(neutFrac05, "", 0.);
0901   INIT_VARIABLE(neutFrac06, "", 0.);
0902   INIT_VARIABLE(neutFrac07, "", 0.);
0903 
0904   INIT_VARIABLE(emFrac01, "", 0.);
0905   INIT_VARIABLE(emFrac02, "", 0.);
0906   INIT_VARIABLE(emFrac03, "", 0.);
0907   INIT_VARIABLE(emFrac04, "", 0.);
0908   INIT_VARIABLE(emFrac05, "", 0.);
0909   INIT_VARIABLE(emFrac06, "", 0.);
0910   INIT_VARIABLE(emFrac07, "", 0.);
0911 
0912   INIT_VARIABLE(beta, "", 0.);
0913   INIT_VARIABLE(betaStar, "", 0.);
0914   INIT_VARIABLE(betaClassic, "", 0.);
0915   INIT_VARIABLE(betaStarClassic, "", 0.);
0916 
0917   INIT_VARIABLE(nvtx, "", 0.);
0918   INIT_VARIABLE(rho, "", 0.);
0919   INIT_VARIABLE(nTrueInt, "", 0.);
0920 
0921   INIT_VARIABLE(jetR, "", 0.);
0922   INIT_VARIABLE(jetRchg, "", 0.);
0923   INIT_VARIABLE(dRMatch, "", 0.);
0924 }
0925 
0926 #undef INIT_VARIABLE