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   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   //Clear everything
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   //Link to the RecoObjects
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     // skip candidates to be ignored in the computation
0066     // of PUPPI's alphas (e.g. electrons and muons if puppiNoLep=True)
0067     if (std::abs(rParticle.id) == 3)
0068       continue;
0069
0070     fPFParticlesForVar.push_back(pCand);
0071     // charged candidates assigned to LV
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 &centre,
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 //In fact takes the median not the average
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     //Calculate the Puppi Algo to use
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     //Get the Puppi Sub Algo (given iteration)
0149     int pAlgo = fPuppiAlgo[pPupId].algoId(iOpt);
0150     bool pCharged = fPuppiAlgo[pPupId].isCharged(iOpt);
0151     double pCone = fPuppiAlgo[pPupId].coneSize(iOpt);
0152     // compute the Puppi metric:
0153     //  - calculate goodVar only for candidates that (1) will not be assigned a predefined weight (e.g 0, 1),
0154     //    or (2) are required for computations inside puppi-algos (see call to PuppiAlgo::add below)
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     // code added by Nhan: now instead for every algorithm give it all the particles
0170     for (int i1 = 0; i1 < fNAlgos; i1++) {
0171       // skip cands outside of algo's etaMaxExtrap, as they would anyway be ignored inside PuppiAlgo::add (see end of the block)
0172       if (not(std::abs(iConstits[i0].eta) < fPuppiAlgo[i1].etaMaxExtrap() and getsDefaultWgtIfApplyCHS))
0173         continue;
0174
0175       auto curVal = pVal;
0176       // recompute goodVar if algo has changed
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 //In fact takes the median not the average
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       //Get the Puppi Sub Algo (given iteration)
0200       int pAlgo = fPuppiAlgo[j0].algoId(iOpt);
0201       bool pCharged = fPuppiAlgo[j0].isCharged(iOpt);
0202       double pCone = fPuppiAlgo[j0].coneSize(iOpt);
0203       //Compute the Puppi Metric
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   //if(lId == -1) std::cerr << "Error : Full fiducial range is not defined " << std::endl;
0230   return lId;
0231 }
0232 double PuppiContainer::getChi2FromdZ(double iDZ) {
0233   //We need to obtain prob of PU + (1-Prob of LV)
0234   // Prob(LV) = Gaus(dZ,sigma) where sigma = 1.5mm  (its really more like 1mm)
0235   //double lProbLV = ROOT::Math::normal_cdf_c(std::abs(iDZ),0.2)*2.; //*2 is to do it double sided
0236   //Take iDZ to be corrected by sigma already
0237   double lProbLV = ROOT::Math::normal_cdf_c(std::abs(iDZ), 1.) * 2.;  //*2 is to do it double sided
0238   double lProbPU = 1 - lProbLV;
0239   if (lProbPU <= 0)
0240     lProbPU = 1e-16;  //Quick Trick to through out infs
0241   if (lProbPU >= 0)
0242     lProbPU = 1 - 1e-16;  //Ditto
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   //Run through all compute mean and RMS
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     //Refresh
0271     pVals.clear();
0272     double pWeight = 1;
0273     //Get the Puppi Id and if ill defined move on
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     // fill the p-values
0284     double pChi2 = 0;
0285     if (fUseExp) {
0286       //Compute an Experimental Puppi Weight with delta Z info (very simple example)
0287       pChi2 = getChi2FromdZ(rParticle.dZ);
0288       //Now make sure Neutrals are not set
0289       if ((std::abs(rParticle.pdgId) == 22) || (std::abs(rParticle.pdgId) == 130))
0290         pChi2 = 0;
0291     }
0292     //Fill and compute the PuppiWeight
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     //Apply the CHS weights
0299     if (rParticle.id == 1 && fApplyCHS)
0300       pWeight = 1;
0301     if (rParticle.id == 2 && fApplyCHS)
0302       pWeight = 0;
0303     //Apply weight of 1 for leptons if puppiNoLep
0304     if (rParticle.id == 3)
0305       pWeight = 1;
0306     //Basic Weight Checks
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     //Basic Cuts
0314     if (pWeight * fPFParticles[i0].pt < fPuppiAlgo[pPupId].neutralPt(fPUProxy) && rParticle.id == 0)
0315       pWeight = 0;  //threshold cut on the neutral Pt
0316     // Protect high pT photons (important for gamma to hadronic recoil balance)
0317     if (fPtMaxPhotons > 0 && rParticle.pdgId == 22 && std::abs(fPFParticles[i0].eta) < fEtaMaxPhotons &&
0318         fPFParticles[i0].pt > fPtMaxPhotons)
0319       pWeight = 1.;
0320     // Protect high pT neutrals
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;  //==> Elminate the low Weight stuff
0326     if (fInvert)
0327       pWeight = 1. - pWeight;
0328     //std::cout << "rParticle.pt = " <<  rParticle.pt << ", rParticle.charge = " << rParticle.charge << ", rParticle.id = " << rParticle.id << ", weight = " << pWeight << std::endl;
0329
0330     fWeights.push_back(pWeight);
0331     fAlphaMed.push_back(fPuppiAlgo[pPupId].median());
0332     fAlphaRMS.push_back(fPuppiAlgo[pPupId].rms());
0333     //Now get rid of the thrown out weights for the particle collection
0334
0335     // leave these lines in, in case want to move eventually to having no 1-to-1 correspondence between puppi and pf cands
0336     // if( std::abs(pWeight) < std::numeric_limits<double>::denorm_min() ) continue; // this line seems not to work like it's supposed to...
0337     // if(std::abs(pWeight) <= 0. ) continue;
0338   }
0339   return fWeights;
0340 }