File indexing completed on 2023-03-17 10:45:24
0001 #include "CommonTools/PileupAlgos/interface/PuppiContainer.h"
0002 #include "DataFormats/Math/interface/deltaR.h"
0003 #include "Math/ProbFunc.h"
0004 #include "TMath.h"
0005 #include <iostream>
0006 #include <cmath>
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "FWCore/Utilities/interface/isFinite.h"
0009
0010 using namespace std;
0011
0012 PuppiContainer::PuppiContainer(const edm::ParameterSet &iConfig) {
0013 fPuppiDiagnostics = iConfig.getParameter<bool>("puppiDiagnostics");
0014 fApplyCHS = iConfig.getParameter<bool>("applyCHS");
0015 fInvert = iConfig.getParameter<bool>("invertPuppi");
0016 fUseExp = iConfig.getParameter<bool>("useExp");
0017 fPuppiWeightCut = iConfig.getParameter<double>("MinPuppiWeight");
0018 fPtMaxPhotons = iConfig.getParameter<double>("PtMaxPhotons");
0019 fEtaMaxPhotons = iConfig.getParameter<double>("EtaMaxPhotons");
0020 fPtMaxNeutrals = iConfig.getParameter<double>("PtMaxNeutrals");
0021 fPtMaxNeutralsStartSlope = iConfig.getParameter<double>("PtMaxNeutralsStartSlope");
0022 std::vector<edm::ParameterSet> lAlgos = iConfig.getParameter<std::vector<edm::ParameterSet> >("algos");
0023 fNAlgos = lAlgos.size();
0024 for (unsigned int i0 = 0; i0 < lAlgos.size(); i0++) {
0025 PuppiAlgo pPuppiConfig(lAlgos[i0]);
0026 fPuppiAlgo.push_back(pPuppiConfig);
0027 }
0028 }
0029
0030 void PuppiContainer::initialize(const std::vector<RecoObj> &iRecoObjects) {
0031
0032 fPFParticles.resize(0);
0033 fPFParticlesForVar.resize(0);
0034 fPFParticlesForVarChargedPV.resize(0);
0035 fWeights.resize(0);
0036 fVals.resize(0);
0037 fRawAlphas.resize(0);
0038 fAlphaMed.resize(0);
0039 fAlphaRMS.resize(0);
0040 fPUProxy = 1.;
0041
0042 fRecoParticles = &iRecoObjects;
0043 fPFParticles.reserve(iRecoObjects.size());
0044 fPFParticlesForVar.reserve(iRecoObjects.size());
0045 fPFParticlesForVarChargedPV.reserve(iRecoObjects.size());
0046 for (auto const &rParticle : *fRecoParticles) {
0047 PuppiCandidate pCand;
0048 pCand.id = rParticle.id;
0049 if (edm::isFinite(rParticle.rapidity)) {
0050 pCand.pt = rParticle.pt;
0051 pCand.eta = rParticle.eta;
0052 pCand.rapidity = rParticle.rapidity;
0053 pCand.phi = rParticle.phi;
0054 pCand.m = rParticle.m;
0055 } else {
0056 pCand.pt = 0.;
0057 pCand.eta = 99.;
0058 pCand.rapidity = 99.;
0059 pCand.phi = 0.;
0060 pCand.m = 0.;
0061 }
0062
0063 fPFParticles.push_back(pCand);
0064
0065
0066
0067 if (std::abs(rParticle.id) == 3)
0068 continue;
0069
0070 fPFParticlesForVar.push_back(pCand);
0071
0072 if (std::abs(rParticle.id) == 1)
0073 fPFParticlesForVarChargedPV.push_back(pCand);
0074 }
0075 }
0076
0077 PuppiContainer::~PuppiContainer() {}
0078
0079 double PuppiContainer::goodVar(PuppiCandidate const &iPart,
0080 std::vector<PuppiCandidate> const &iParts,
0081 int iOpt,
0082 const double iRCone) {
0083 return var_within_R(iOpt, iParts, iPart, iRCone);
0084 }
0085
0086 double PuppiContainer::var_within_R(int iId,
0087 const vector<PuppiCandidate> &particles,
0088 const PuppiCandidate ¢re,
0089 const double R) {
0090 if (iId == -1)
0091 return 1.;
0092
0093 double const r2 = R * R;
0094 double var = 0.;
0095
0096 for (auto const &cand : particles) {
0097 if (std::abs(cand.rapidity - centre.rapidity) < R) {
0098 auto const dr2y = reco::deltaR2(cand.rapidity, cand.phi, centre.rapidity, centre.phi);
0099 if (dr2y < r2) {
0100 auto const dr2 = reco::deltaR2(cand.eta, cand.phi, centre.eta, centre.phi);
0101 if (dr2 < 0.0001)
0102 continue;
0103 auto const pt = cand.pt;
0104 switch (iId) {
0105 case 5:
0106 var += (pt * pt / dr2);
0107 break;
0108 case 4:
0109 var += pt;
0110 break;
0111 case 3:
0112 var += (1. / dr2);
0113 break;
0114 case 2:
0115 var += (1. / dr2);
0116 break;
0117 case 1:
0118 var += pt;
0119 break;
0120 case 0:
0121 var += (pt / dr2);
0122 break;
0123 }
0124 }
0125 }
0126 }
0127
0128 if ((var != 0.) and ((iId == 0) or (iId == 3) or (iId == 5)))
0129 var = log(var);
0130 else if (iId == 1)
0131 var += centre.pt;
0132
0133 return var;
0134 }
0135
0136
0137 void PuppiContainer::getRMSAvg(int iOpt,
0138 std::vector<PuppiCandidate> const &iConstits,
0139 std::vector<PuppiCandidate> const &iParticles,
0140 std::vector<PuppiCandidate> const &iChargedParticles) {
0141 for (unsigned int i0 = 0; i0 < iConstits.size(); i0++) {
0142
0143 int pPupId = getPuppiId(iConstits[i0].pt, iConstits[i0].eta);
0144 if (pPupId == -1 || fPuppiAlgo[pPupId].numAlgos() <= iOpt) {
0145 fVals.push_back(-1);
0146 continue;
0147 }
0148
0149 int pAlgo = fPuppiAlgo[pPupId].algoId(iOpt);
0150 bool pCharged = fPuppiAlgo[pPupId].isCharged(iOpt);
0151 double pCone = fPuppiAlgo[pPupId].coneSize(iOpt);
0152
0153
0154
0155 double pVal = -1;
0156 bool const getsDefaultWgtIfApplyCHS = iConstits[i0].id == 1 or iConstits[i0].id == 2;
0157 if (not((fApplyCHS and getsDefaultWgtIfApplyCHS) or iConstits[i0].id == 3) or
0158 (std::abs(iConstits[i0].eta) < fPuppiAlgo[pPupId].etaMaxExtrap() and getsDefaultWgtIfApplyCHS)) {
0159 pVal = goodVar(iConstits[i0], pCharged ? iChargedParticles : iParticles, pAlgo, pCone);
0160 }
0161 fVals.push_back(pVal);
0162
0163 if (!edm::isFinite(pVal)) {
0164 LogDebug("NotFound") << "====> Value is Nan " << pVal << " == " << iConstits[i0].pt << " -- " << iConstits[i0].eta
0165 << endl;
0166 continue;
0167 }
0168
0169
0170 for (int i1 = 0; i1 < fNAlgos; i1++) {
0171
0172 if (not(std::abs(iConstits[i0].eta) < fPuppiAlgo[i1].etaMaxExtrap() and getsDefaultWgtIfApplyCHS))
0173 continue;
0174
0175 auto curVal = pVal;
0176
0177 if (i1 != pPupId) {
0178 pAlgo = fPuppiAlgo[i1].algoId(iOpt);
0179 pCharged = fPuppiAlgo[i1].isCharged(iOpt);
0180 pCone = fPuppiAlgo[i1].coneSize(iOpt);
0181 curVal = goodVar(iConstits[i0], pCharged ? iChargedParticles : iParticles, pAlgo, pCone);
0182 }
0183
0184 fPuppiAlgo[i1].add(iConstits[i0], curVal, iOpt);
0185 }
0186 }
0187
0188 for (int i0 = 0; i0 < fNAlgos; i0++)
0189 fPuppiAlgo[i0].computeMedRMS(iOpt);
0190 }
0191
0192
0193 void PuppiContainer::getRawAlphas(int iOpt,
0194 std::vector<PuppiCandidate> const &iConstits,
0195 std::vector<PuppiCandidate> const &iParticles,
0196 std::vector<PuppiCandidate> const &iChargedParticles) {
0197 for (int j0 = 0; j0 < fNAlgos; j0++) {
0198 for (unsigned int i0 = 0; i0 < iConstits.size(); i0++) {
0199
0200 int pAlgo = fPuppiAlgo[j0].algoId(iOpt);
0201 bool pCharged = fPuppiAlgo[j0].isCharged(iOpt);
0202 double pCone = fPuppiAlgo[j0].coneSize(iOpt);
0203
0204 double const pVal = goodVar(iConstits[i0], pCharged ? iChargedParticles : iParticles, pAlgo, pCone);
0205 fRawAlphas.push_back(pVal);
0206 if (!edm::isFinite(pVal)) {
0207 LogDebug("NotFound") << "====> Value is Nan " << pVal << " == " << iConstits[i0].pt << " -- "
0208 << iConstits[i0].eta << endl;
0209 continue;
0210 }
0211 }
0212 }
0213 }
0214
0215 int PuppiContainer::getPuppiId(float iPt, float iEta) {
0216 int lId = -1;
0217 for (int i0 = 0; i0 < fNAlgos; i0++) {
0218 int nEtaBinsPerAlgo = fPuppiAlgo[i0].etaBins();
0219 for (int i1 = 0; i1 < nEtaBinsPerAlgo; i1++) {
0220 if ((std::abs(iEta) >= fPuppiAlgo[i0].etaMin(i1)) && (std::abs(iEta) < fPuppiAlgo[i0].etaMax(i1))) {
0221 fPuppiAlgo[i0].fixAlgoEtaBin(i1);
0222 if (iPt > fPuppiAlgo[i0].ptMin()) {
0223 lId = i0;
0224 break;
0225 }
0226 }
0227 }
0228 }
0229
0230 return lId;
0231 }
0232 double PuppiContainer::getChi2FromdZ(double iDZ) {
0233
0234
0235
0236
0237 double lProbLV = ROOT::Math::normal_cdf_c(std::abs(iDZ), 1.) * 2.;
0238 double lProbPU = 1 - lProbLV;
0239 if (lProbPU <= 0)
0240 lProbPU = 1e-16;
0241 if (lProbPU >= 0)
0242 lProbPU = 1 - 1e-16;
0243 double lChi2PU = TMath::ChisquareQuantile(lProbPU, 1);
0244 lChi2PU *= lChi2PU;
0245 return lChi2PU;
0246 }
0247 std::vector<double> const &PuppiContainer::puppiWeights() {
0248 int lNParticles = fRecoParticles->size();
0249
0250 fWeights.clear();
0251 fWeights.reserve(lNParticles);
0252 fVals.clear();
0253 fVals.reserve(lNParticles);
0254 for (int i0 = 0; i0 < fNAlgos; i0++)
0255 fPuppiAlgo[i0].reset();
0256
0257 int lNMaxAlgo = 1;
0258 for (int i0 = 0; i0 < fNAlgos; i0++)
0259 lNMaxAlgo = std::max(fPuppiAlgo[i0].numAlgos(), lNMaxAlgo);
0260
0261 for (int i0 = 0; i0 < lNMaxAlgo; i0++) {
0262 getRMSAvg(i0, fPFParticles, fPFParticlesForVar, fPFParticlesForVarChargedPV);
0263 }
0264 if (fPuppiDiagnostics)
0265 getRawAlphas(0, fPFParticles, fPFParticlesForVar, fPFParticlesForVarChargedPV);
0266
0267 std::vector<double> pVals;
0268 pVals.reserve(lNParticles);
0269 for (int i0 = 0; i0 < lNParticles; i0++) {
0270
0271 pVals.clear();
0272 double pWeight = 1;
0273
0274 const auto &rParticle = (*fRecoParticles)[i0];
0275 int pPupId = getPuppiId(rParticle.pt, rParticle.eta);
0276 if (pPupId == -1) {
0277 fWeights.push_back(0);
0278 fAlphaMed.push_back(-10);
0279 fAlphaRMS.push_back(-10);
0280 continue;
0281 }
0282
0283
0284 double pChi2 = 0;
0285 if (fUseExp) {
0286
0287 pChi2 = getChi2FromdZ(rParticle.dZ);
0288
0289 if ((std::abs(rParticle.pdgId) == 22) || (std::abs(rParticle.pdgId) == 130))
0290 pChi2 = 0;
0291 }
0292
0293 int lNAlgos = fPuppiAlgo[pPupId].numAlgos();
0294 for (int i1 = 0; i1 < lNAlgos; i1++)
0295 pVals.push_back(fVals[lNParticles * i1 + i0]);
0296
0297 pWeight = fPuppiAlgo[pPupId].compute(pVals, pChi2);
0298
0299 if (rParticle.id == 1 && fApplyCHS)
0300 pWeight = 1;
0301 if (rParticle.id == 2 && fApplyCHS)
0302 pWeight = 0;
0303
0304 if (rParticle.id == 3)
0305 pWeight = 1;
0306
0307 if (!edm::isFinite(pWeight)) {
0308 pWeight = 0.0;
0309 LogDebug("PuppiWeightError") << "====> Weight is nan : " << pWeight << " : pt " << rParticle.pt
0310 << " -- eta : " << rParticle.eta << " -- Value" << fVals[i0]
0311 << " -- id : " << rParticle.id << " -- NAlgos: " << lNAlgos << std::endl;
0312 }
0313
0314 if (pWeight * fPFParticles[i0].pt < fPuppiAlgo[pPupId].neutralPt(fPUProxy) && rParticle.id == 0)
0315 pWeight = 0;
0316
0317 if (fPtMaxPhotons > 0 && rParticle.pdgId == 22 && std::abs(fPFParticles[i0].eta) < fEtaMaxPhotons &&
0318 fPFParticles[i0].pt > fPtMaxPhotons)
0319 pWeight = 1.;
0320
0321 else if ((fPtMaxNeutrals > 0) && (rParticle.id == 0))
0322 pWeight = std::clamp(
0323 (fPFParticles[i0].pt - fPtMaxNeutralsStartSlope) / (fPtMaxNeutrals - fPtMaxNeutralsStartSlope), pWeight, 1.);
0324 if (pWeight < fPuppiWeightCut)
0325 pWeight = 0;
0326 if (fInvert)
0327 pWeight = 1. - pWeight;
0328
0329
0330 fWeights.push_back(pWeight);
0331 fAlphaMed.push_back(fPuppiAlgo[pPupId].median());
0332 fAlphaRMS.push_back(fPuppiAlgo[pPupId].rms());
0333
0334
0335
0336
0337
0338 }
0339 return fWeights;
0340 }