Phase2L1ParticleFlow/src/LinearizedPuppiAlgo.cc

0001 #include "L1Trigger/Phase2L1ParticleFlow/interface/LinearizedPuppiAlgo.h"
0002 #include "DataFormats/L1TParticleFlow/interface/PFCandidate.h"
0003 #include "FWCore/Utilities/interface/Exception.h"
0004 #include "DataFormats/Math/interface/deltaR.h"
0005 #include "L1Trigger/Phase2L1ParticleFlow/interface/dbgPrintf.h"
0006
0007 #include "Math/ProbFunc.h"
0008
0009 namespace {
0010   std::vector<float> vd2vf(const std::vector<double> &vd) {
0011     std::vector<float> ret;
0012     ret.insert(ret.end(), vd.begin(), vd.end());
0013     return ret;
0014   }
0015 }  // namespace
0016
0017 using namespace l1tpf_impl;
0018
0019 LinearizedPuppiAlgo::LinearizedPuppiAlgo(const edm::ParameterSet &iConfig)
0020     : PuppiAlgo(iConfig),
0021       puppiPriors_(vd2vf(iConfig.getParameter<std::vector<double>>("puppiPriors"))),
0022       puppiPriorsPhotons_(vd2vf(iConfig.getParameter<std::vector<double>>("puppiPriorsPhotons"))),
0023       puppiPtSlopes_(vd2vf(iConfig.getParameter<std::vector<double>>("puppiPtSlopes"))),
0024       puppiPtSlopesPhotons_(vd2vf(iConfig.getParameter<std::vector<double>>("puppiPtSlopesPhotons"))),
0025       puppiPtZeros_(vd2vf(iConfig.getParameter<std::vector<double>>("puppiPtZeros"))),
0026       puppiPtZerosPhotons_(vd2vf(iConfig.getParameter<std::vector<double>>("puppiPtZerosPhotons"))),
0027       puppiAlphaSlopes_(vd2vf(iConfig.getParameter<std::vector<double>>("puppiAlphaSlopes"))),
0028       puppiAlphaSlopesPhotons_(vd2vf(iConfig.getParameter<std::vector<double>>("puppiAlphaSlopesPhotons"))),
0029       puppiAlphaZeros_(vd2vf(iConfig.getParameter<std::vector<double>>("puppiAlphaZeros"))),
0030       puppiAlphaZerosPhotons_(vd2vf(iConfig.getParameter<std::vector<double>>("puppiAlphaZerosPhotons"))),
0031       puppiAlphaCrops_(vd2vf(iConfig.getParameter<std::vector<double>>("puppiAlphaCrops"))),
0032       puppiAlphaCropsPhotons_(vd2vf(iConfig.getParameter<std::vector<double>>("puppiAlphaCropsPhotons"))) {
0033   if (puppiPriors_.size() != puppiEtaCuts_.size())
0034     throw cms::Exception("Configuration", "Mismatched lenght for puppiPriors\n");
0035   if (puppiPtSlopes_.size() != puppiEtaCuts_.size())
0036     throw cms::Exception("Configuration", "Mismatched lenght for puppiPtSlopes\n");
0037   if (puppiPtZeros_.size() != puppiEtaCuts_.size())
0038     throw cms::Exception("Configuration", "Mismatched lenght for puppiPtZeros\n");
0039   if (puppiAlphaSlopes_.size() != puppiEtaCuts_.size())
0040     throw cms::Exception("Configuration", "Mismatched lenght for puppiAlphaSlopes\n");
0041   if (puppiAlphaZeros_.size() != puppiEtaCuts_.size())
0042     throw cms::Exception("Configuration", "Mismatched lenght for puppiAlphaZeros\n");
0043   if (puppiAlphaCrops_.size() != puppiEtaCuts_.size())
0044     throw cms::Exception("Configuration", "Mismatched lenght for puppiAlphaCrops\n");
0045   if (puppiPriorsPhotons_.size() != puppiEtaCuts_.size())
0046     throw cms::Exception("Configuration", "Mismatched lenght for puppiPriorsPhotons\n");
0047   if (puppiPtSlopesPhotons_.size() != puppiEtaCuts_.size())
0048     throw cms::Exception("Configuration", "Mismatched lenght for puppiPtSlopesPhotons\n");
0049   if (puppiPtZerosPhotons_.size() != puppiEtaCuts_.size())
0050     throw cms::Exception("Configuration", "Mismatched lenght for puppiPtZerosPhotons\n");
0051   if (puppiAlphaSlopesPhotons_.size() != puppiEtaCuts_.size())
0052     throw cms::Exception("Configuration", "Mismatched lenght for puppiAlphaSlopesPhotons\n");
0053   if (puppiAlphaZerosPhotons_.size() != puppiEtaCuts_.size())
0054     throw cms::Exception("Configuration", "Mismatched lenght for puppiAlphaZerosPhotons\n");
0055   if (puppiAlphaCropsPhotons_.size() != puppiEtaCuts_.size())
0056     throw cms::Exception("Configuration", "Mismatched lenght for puppiAlphaCropsPhotons\n");
0057 }
0058
0059 LinearizedPuppiAlgo::~LinearizedPuppiAlgo() {}
0060
0061 const std::vector<std::string> &LinearizedPuppiAlgo::puGlobalNames() const {
0062   static const std::vector<std::string> names_{};
0063   return names_;
0064 }
0065 void LinearizedPuppiAlgo::doPUGlobals(const std::vector<Region> &rs,
0066                                       float z0,
0067                                       float npu,
0068                                       std::vector<float> &globals) const {
0069   globals.clear();
0070 }
0071 void LinearizedPuppiAlgo::runNeutralsPU(Region &r, float z0, float npu, const std::vector<float> &globals) const {
0072   std::vector<float> alphaC, alphaF;
0073   PuppiAlgo::computePuppiAlphas(r, alphaC, alphaF);
0074   computePuppiWeights(r, npu, alphaC, alphaF);
0075   PuppiAlgo::fillPuppi(r);
0076 }
0077
0078 void LinearizedPuppiAlgo::computePuppiWeights(Region &r,
0079                                               float npu,
0080                                               const std::vector<float> &alphaC,
0081                                               const std::vector<float> &alphaF) const {
0082   if (debug_ && npu > 0)
0083     dbgPrintf("LinPup\t npu estimate %7.2f --> log(npu/200) = %+6.2f \n", npu, std::log(npu / 200.f));
0084   for (unsigned int ip = 0, np = r.pf.size(); ip < np; ++ip) {
0085     PFParticle &p = r.pf[ip];
0086     // charged
0087     if (p.hwId == l1t::PFCandidate::ChargedHadron || p.hwId == l1t::PFCandidate::Electron ||
0088         p.hwId == l1t::PFCandidate::Muon) {
0089       p.setPuppiW(p.chargedPV || p.hwId == l1t::PFCandidate::Muon ? 1.0 : 0);
0090       if (debug_ == 2)
0091         dbgPrintf(
0092             "LinPup\t charged id %1d pt %7.2f eta %+5.2f phi %+5.2f   fromPV %1d                                       "
0093             "                        --> puppi weight %.3f   puppi pt %7.2f \n",
0094             p.hwId,
0095             p.floatPt(),
0096             p.floatEta(),
0097             p.floatPhi(),
0098             p.chargedPV,
0099             p.floatPuppiW(),
0100             p.floatPt() * p.floatPuppiW());
0101       continue;
0102     }
0103     // neutral
0104     float absEta = r.relativeCoordinates ? r.globalAbsEta(p.floatEta()) : std::abs(p.floatEta());
0105     bool central = absEta < etaCharged_;  // FIXME could make a better integer implementation
0106     bool photon = (p.hwId == l1t::PFCandidate::Photon);
0107     // get alpha
0108     float alpha = central ? alphaC[ip] : alphaF[ip];
0109     alpha = (alpha > 0 ? std::log(alpha) : 0);
0110     // get eta bin
0111     unsigned int ietaBin = 0, lastBin = puppiEtaCuts_.size() - 1;
0112     while (ietaBin < lastBin && absEta > puppiEtaCuts_[ietaBin]) {
0113       ietaBin++;
0114     }
0115     float alphaZero = (photon ? puppiAlphaZerosPhotons_ : puppiAlphaZeros_)[ietaBin];
0116     float alphaSlope = (photon ? puppiAlphaSlopesPhotons_ : puppiAlphaSlopes_)[ietaBin];
0117     float alphaCrop = (photon ? puppiAlphaCropsPhotons_ : puppiAlphaCrops_)[ietaBin];
0118     float x2a = std::clamp(alphaSlope * (alpha - alphaZero), -alphaCrop, alphaCrop);
0119     // weight by pT
0120     float ptZero = (photon ? puppiPtZerosPhotons_ : puppiPtZeros_)[ietaBin];
0121     float ptSlope = (photon ? puppiPtSlopesPhotons_ : puppiPtSlopes_)[ietaBin];
0122     float x2pt = ptSlope * (p.floatPt() - ptZero);
0123     // weight by prior
0124     float prior = (photon ? puppiPriorsPhotons_ : puppiPriors_)[ietaBin];
0125     float x2prior = (npu > 0 ? std::log(npu / 200.f) : 0) + prior;
0126     // total
0127     float x2 = x2a + x2pt - x2prior;
0128     p.setPuppiW(1.0 / (1.0 + std::exp(-x2)));
0129     if (debug_ == 1 || debug_ == 2 || debug_ == int(10 + ietaBin))
0130       dbgPrintf(
0131           "LinPup\t neutral id %1d pt %7.2f eta %+5.2f phi %+5.2f   alpha %+6.2f   x2a %+5.2f   x2pt %+6.2f   x2prior "
0132           "%+6.2f -->  x2 %+6.2f --> puppi weight %.3f   puppi pt %7.2f \n",
0133           p.hwId,
0134           p.floatPt(),
0135           p.floatEta(),
0136           p.floatPhi(),
0137           alpha,
0138           x2a,
0139           x2pt,
0140           -x2prior,
0141           x2,
0142           p.floatPuppiW(),
0143           p.floatPt() * p.floatPuppiW());
0144   }
0145 }