PileupAlgos/src/PuppiContainer.cc

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   fPuppiAlgo.reserve(lAlgos.size());
0024   for (auto const &algos : lAlgos) {
0025     fPuppiAlgo.emplace_back(algos);
0026   }
0027 }
0028
0029 void PuppiContainer::initialize(const std::vector<RecoObj> &iRecoObjects,
0030                                 std::vector<PuppiCandidate> &pfParticles,
0031                                 std::vector<PuppiCandidate> &pfParticlesForVar,
0032                                 std::vector<PuppiCandidate> &pfParticlesForVarChargedPV) const {
0033   pfParticles.reserve(iRecoObjects.size());
0034   pfParticlesForVar.reserve(iRecoObjects.size());
0035   pfParticlesForVarChargedPV.reserve(iRecoObjects.size());
0036   for (auto const &rParticle : iRecoObjects) {
0037     PuppiCandidate pCand;
0038     pCand.id = rParticle.id;
0039     if (edm::isFinite(rParticle.rapidity)) {
0040       pCand.pt = rParticle.pt;
0041       pCand.eta = rParticle.eta;
0042       pCand.rapidity = rParticle.rapidity;
0043       pCand.phi = rParticle.phi;
0044       pCand.m = rParticle.m;
0045     } else {
0046       pCand.pt = 0.;
0047       pCand.eta = 99.;
0048       pCand.rapidity = 99.;
0049       pCand.phi = 0.;
0050       pCand.m = 0.;
0051     }
0052
0053     pfParticles.push_back(pCand);
0054
0055     // skip candidates to be ignored in the computation
0056     // of PUPPI's alphas (e.g. electrons and muons if puppiNoLep=True)
0057     if (std::abs(rParticle.id) == 3)
0058       continue;
0059
0060     pfParticlesForVar.push_back(pCand);
0061     // charged candidates assigned to LV
0062     if (std::abs(rParticle.id) == 1)
0063       pfParticlesForVarChargedPV.push_back(pCand);
0064   }
0065 }
0066
0067 double PuppiContainer::goodVar(PuppiCandidate const &iPart,
0068                                std::vector<PuppiCandidate> const &iParts,
0069                                int iOpt,
0070                                const double iRCone) const {
0071   return var_within_R(iOpt, iParts, iPart, iRCone);
0072 }
0073
0074 double PuppiContainer::var_within_R(int iId,
0075                                     const vector<PuppiCandidate> &particles,
0076                                     const PuppiCandidate &centre,
0077                                     const double R) const {
0078   if (iId == -1)
0079     return 1.;
0080
0081   double const r2 = R * R;
0082   double var = 0.;
0083
0084   for (auto const &cand : particles) {
0085     if (std::abs(cand.rapidity - centre.rapidity) < R) {
0086       auto const dr2y = reco::deltaR2(cand.rapidity, cand.phi, centre.rapidity, centre.phi);
0087       if (dr2y < r2) {
0088         auto const dr2 = reco::deltaR2(cand.eta, cand.phi, centre.eta, centre.phi);
0089         if (dr2 < 0.0001)
0090           continue;
0091         auto const pt = cand.pt;
0092         switch (iId) {
0093           case 5:
0094             var += (pt * pt / dr2);
0095             break;
0096           case 4:
0097             var += pt;
0098             break;
0099           case 3:
0100             var += (1. / dr2);
0101             break;
0102           case 2:
0103             var += (1. / dr2);
0104             break;
0105           case 1:
0106             var += pt;
0107             break;
0108           case 0:
0109             var += (pt / dr2);
0110             break;
0111         }
0112       }
0113     }
0114   }
0115
0116   if ((var != 0.) and ((iId == 0) or (iId == 3) or (iId == 5)))
0117     var = log(var);
0118   else if (iId == 1)
0119     var += centre.pt;
0120
0121   return var;
0122 }
0123
0124 //In fact takes the median not the average
0125 void PuppiContainer::getRMSAvg(int iOpt,
0126                                std::vector<PuppiCandidate> const &iConstits,
0127                                std::vector<PuppiCandidate> const &iParticles,
0128                                std::vector<PuppiCandidate> const &iChargedParticles,
0129                                std::vector<double> &oVals) {
0130   for (unsigned int i0 = 0; i0 < iConstits.size(); i0++) {
0131     //Calculate the Puppi Algo to use
0132     int pPupId = getPuppiId(iConstits[i0].pt, iConstits[i0].eta);
0133     if (pPupId == -1 || fPuppiAlgo[pPupId].numAlgos() <= iOpt) {
0134       oVals.push_back(-1);
0135       continue;
0136     }
0137     //Get the Puppi Sub Algo (given iteration)
0138     int pAlgo = fPuppiAlgo[pPupId].algoId(iOpt);
0139     bool pCharged = fPuppiAlgo[pPupId].isCharged(iOpt);
0140     double pCone = fPuppiAlgo[pPupId].coneSize(iOpt);
0141     // compute the Puppi metric:
0142     //  - calculate goodVar only for candidates that (1) will not be assigned a predefined weight (e.g 0, 1),
0143     //    or (2) are required for computations inside puppi-algos (see call to PuppiAlgo::add below)
0144     double pVal = -1;
0145     bool const getsDefaultWgtIfApplyCHS = iConstits[i0].id == 1 or iConstits[i0].id == 2;
0146     if (not((fApplyCHS and getsDefaultWgtIfApplyCHS) or iConstits[i0].id == 3) or
0147         (std::abs(iConstits[i0].eta) < fPuppiAlgo[pPupId].etaMaxExtrap() and getsDefaultWgtIfApplyCHS)) {
0148       pVal = goodVar(iConstits[i0], pCharged ? iChargedParticles : iParticles, pAlgo, pCone);
0149     }
0150     oVals.push_back(pVal);
0151
0152     if (!edm::isFinite(pVal)) {
0153       LogDebug("NotFound") << "====> Value is Nan " << pVal << " == " << iConstits[i0].pt << " -- " << iConstits[i0].eta
0154                            << endl;
0155       continue;
0156     }
0157
0158     // code added by Nhan: now instead for every algorithm give it all the particles
0159     int count = 0;
0160     for (auto &algo : fPuppiAlgo) {
0161       int index = count++;
0162       // skip cands outside of algo's etaMaxExtrap, as they would anyway be ignored inside PuppiAlgo::add (see end of the block)
0163       if (not(std::abs(iConstits[i0].eta) < algo.etaMaxExtrap() and getsDefaultWgtIfApplyCHS))
0164         continue;
0165
0166       auto curVal = pVal;
0167       // recompute goodVar if algo has changed
0168       if (index != pPupId) {
0169         pAlgo = algo.algoId(iOpt);
0170         pCharged = algo.isCharged(iOpt);
0171         pCone = algo.coneSize(iOpt);
0172         curVal = goodVar(iConstits[i0], pCharged ? iChargedParticles : iParticles, pAlgo, pCone);
0173       }
0174
0175       algo.add(iConstits[i0], curVal, iOpt);
0176     }
0177   }
0178
0179   for (auto &algo : fPuppiAlgo)
0180     algo.computeMedRMS(iOpt);
0181 }
0182
0183 //In fact takes the median not the average
0184 std::vector<double> PuppiContainer::getRawAlphas(int iOpt,
0185                                                  std::vector<PuppiCandidate> const &iConstits,
0186                                                  std::vector<PuppiCandidate> const &iParticles,
0187                                                  std::vector<PuppiCandidate> const &iChargedParticles) const {
0188   std::vector<double> oRawAlphas;
0189   oRawAlphas.reserve(fPuppiAlgo.size() * iConstits.size());
0190   for (auto &algo : fPuppiAlgo) {
0191     for (auto const &constit : iConstits) {
0192       //Get the Puppi Sub Algo (given iteration)
0193       int pAlgo = algo.algoId(iOpt);
0194       bool pCharged = algo.isCharged(iOpt);
0195       double pCone = algo.coneSize(iOpt);
0196       //Compute the Puppi Metric
0197       double const pVal = goodVar(constit, pCharged ? iChargedParticles : iParticles, pAlgo, pCone);
0198       oRawAlphas.push_back(pVal);
0199       if (!edm::isFinite(pVal)) {
0200         LogDebug("NotFound") << "====> Value is Nan " << pVal << " == " << constit.pt << " -- " << constit.eta << endl;
0201         continue;
0202       }
0203     }
0204   }
0205   return oRawAlphas;
0206 }
0207
0208 int PuppiContainer::getPuppiId(float iPt, float iEta) {
0209   int lId = -1;
0210   int count = 0;
0211   for (auto &algo : fPuppiAlgo) {
0212     int index = count++;
0213     int nEtaBinsPerAlgo = algo.etaBins();
0214     for (int i1 = 0; i1 < nEtaBinsPerAlgo; i1++) {
0215       if ((std::abs(iEta) >= algo.etaMin(i1)) && (std::abs(iEta) < algo.etaMax(i1))) {
0216         algo.fixAlgoEtaBin(i1);
0217         if (iPt > algo.ptMin()) {
0218           lId = index;
0219           break;
0220         }
0221       }
0222     }
0223   }
0224   //if(lId == -1) std::cerr << "Error : Full fiducial range is not defined " << std::endl;
0225   return lId;
0226 }
0227 double PuppiContainer::getChi2FromdZ(double iDZ) const {
0228   //We need to obtain prob of PU + (1-Prob of LV)
0229   // Prob(LV) = Gaus(dZ,sigma) where sigma = 1.5mm  (its really more like 1mm)
0230   //double lProbLV = ROOT::Math::normal_cdf_c(std::abs(iDZ),0.2)*2.; //*2 is to do it double sided
0231   //Take iDZ to be corrected by sigma already
0232   double lProbLV = ROOT::Math::normal_cdf_c(std::abs(iDZ), 1.) * 2.;  //*2 is to do it double sided
0233   double lProbPU = 1 - lProbLV;
0234   if (lProbPU <= 0)
0235     lProbPU = 1e-16;  //Quick Trick to through out infs
0236   if (lProbPU >= 0)
0237     lProbPU = 1 - 1e-16;  //Ditto
0238   double lChi2PU = TMath::ChisquareQuantile(lProbPU, 1);
0239   lChi2PU *= lChi2PU;
0240   return lChi2PU;
0241 }
0242 PuppiContainer::Weights PuppiContainer::calculatePuppiWeights(const std::vector<RecoObj> &iRecoObjects,
0243                                                               double iPUProxy) {
0244   std::vector<PuppiCandidate> pfParticles;
0245   std::vector<PuppiCandidate> pfParticlesForVar;
0246   std::vector<PuppiCandidate> pfParticlesForVarChargedPV;
0247
0248   initialize(iRecoObjects, pfParticles, pfParticlesForVar, pfParticlesForVarChargedPV);
0249
0250   int lNParticles = iRecoObjects.size();
0251
0252   Weights returnValue;
0253   returnValue.weights.reserve(lNParticles);
0254   returnValue.puppiAlphas.reserve(lNParticles);
0255
0256   //guarantee all algos are rest before leaving this function
0257   auto doReset = [this](void *) {
0258     for (auto &algo : fPuppiAlgo)
0259       algo.reset();
0260   };
0261   std::unique_ptr<decltype(fPuppiAlgo), decltype(doReset)> guard(&fPuppiAlgo, doReset);
0262
0263   int lNMaxAlgo = 1;
0264   for (auto &algo : fPuppiAlgo)
0265     lNMaxAlgo = std::max(algo.numAlgos(), lNMaxAlgo);
0266   //Run through all compute mean and RMS
0267   for (int i0 = 0; i0 < lNMaxAlgo; i0++) {
0268     getRMSAvg(i0, pfParticles, pfParticlesForVar, pfParticlesForVarChargedPV, returnValue.puppiAlphas);
0269   }
0270   if (fPuppiDiagnostics)
0271     returnValue.puppiRawAlphas = getRawAlphas(0, pfParticles, pfParticlesForVar, pfParticlesForVarChargedPV);
0272
0273   std::vector<double> pVals;
0274   pVals.reserve(lNParticles);
0275   for (int i0 = 0; i0 < lNParticles; i0++) {
0276     //Refresh
0277     pVals.clear();
0278     //Get the Puppi Id and if ill defined move on
0279     const auto &rParticle = iRecoObjects[i0];
0280     int pPupId = getPuppiId(rParticle.pt, rParticle.eta);
0281     if (pPupId == -1) {
0282       returnValue.weights.push_back(0);
0283       returnValue.puppiAlphasMed.push_back(-10);
0284       returnValue.puppiAlphasRMS.push_back(-10);
0285       continue;
0286     }
0287
0288     // fill the p-values
0289     double pChi2 = 0;
0290     if (fUseExp) {
0291       //Compute an Experimental Puppi Weight with delta Z info (very simple example)
0292       pChi2 = getChi2FromdZ(rParticle.dZ);
0293       //Now make sure Neutrals are not set
0294       if ((std::abs(rParticle.pdgId) == 22) || (std::abs(rParticle.pdgId) == 130))
0295         pChi2 = 0;
0296     }
0297     //Fill and compute the PuppiWeight
0298     int lNAlgos = fPuppiAlgo[pPupId].numAlgos();
0299     for (int i1 = 0; i1 < lNAlgos; i1++)
0300       pVals.push_back(returnValue.puppiAlphas[lNParticles * i1 + i0]);
0301
0302     double pWeight = fPuppiAlgo[pPupId].compute(pVals, pChi2);
0303     //Apply the CHS weights
0304     if (rParticle.id == 1 && fApplyCHS)
0305       pWeight = 1;
0306     if (rParticle.id == 2 && fApplyCHS)
0307       pWeight = 0;
0308     //Apply weight of 1 for leptons if puppiNoLep
0309     if (rParticle.id == 3)
0310       pWeight = 1;
0311     //Basic Weight Checks
0312     if (!edm::isFinite(pWeight)) {
0313       pWeight = 0.0;
0314       LogDebug("PuppiWeightError") << "====> Weight is nan : " << pWeight << " : pt " << rParticle.pt
0315                                    << " -- eta : " << rParticle.eta << " -- Value" << returnValue.puppiAlphas[i0]
0316                                    << " -- id :  " << rParticle.id << " --  NAlgos: " << lNAlgos << std::endl;
0317     }
0318     //Basic Cuts
0319     if (pWeight * pfParticles[i0].pt < fPuppiAlgo[pPupId].neutralPt(iPUProxy) && rParticle.id == 0)
0320       pWeight = 0;  //threshold cut on the neutral Pt
0321     // Protect high pT photons (important for gamma to hadronic recoil balance)
0322     if (fPtMaxPhotons > 0 && rParticle.pdgId == 22 && std::abs(pfParticles[i0].eta) < fEtaMaxPhotons &&
0323         pfParticles[i0].pt > fPtMaxPhotons)
0324       pWeight = 1.;
0325     // Protect high pT neutrals
0326     else if ((fPtMaxNeutrals > 0) && (rParticle.id == 0))
0327       pWeight = std::clamp(
0328           (pfParticles[i0].pt - fPtMaxNeutralsStartSlope) / (fPtMaxNeutrals - fPtMaxNeutralsStartSlope), pWeight, 1.);
0329     if (pWeight < fPuppiWeightCut)
0330       pWeight = 0;  //==> Elminate the low Weight stuff
0331     if (fInvert)
0332       pWeight = 1. - pWeight;
0333     //std::cout << "rParticle.pt = " <<  rParticle.pt << ", rParticle.charge = " << rParticle.charge << ", rParticle.id = " << rParticle.id << ", weight = " << pWeight << std::endl;
0334
0335     returnValue.weights.push_back(pWeight);
0336     returnValue.puppiAlphasMed.push_back(fPuppiAlgo[pPupId].median());
0337     returnValue.puppiAlphasRMS.push_back(fPuppiAlgo[pPupId].rms());
0338     //Now get rid of the thrown out returnValue.weights for the particle collection
0339
0340     // leave these lines in, in case want to move eventually to having no 1-to-1 correspondence between puppi and pf cands
0341     // if( std::abs(pWeight) < std::numeric_limits<double>::denorm_min() ) continue; // this line seems not to work like it's supposed to...
0342     // if(std::abs(pWeight) <= 0. ) continue;
0343   }
0344   return returnValue;
0345 }