File indexing completed on 2022-06-10 01:54:01
0001 #include "L1Trigger/Phase2L1ParticleFlow/interface/PuppiAlgo.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 }
0016
0017 using namespace l1tpf_impl;
0018
0019 PuppiAlgo::PuppiAlgo(const edm::ParameterSet &iConfig)
0020 : PUAlgoBase(iConfig),
0021 puppiDr_(iConfig.getParameter<double>("puppiDr")),
0022 puppiDrMin_(iConfig.getParameter<double>("puppiDrMin")),
0023 puppiPtMax_(iConfig.getParameter<double>("puppiPtMax")),
0024 puppiEtaCuts_(vd2vf(iConfig.getParameter<std::vector<double>>("puppiEtaCuts"))),
0025 puppiPtCuts_(vd2vf(iConfig.getParameter<std::vector<double>>("puppiPtCuts"))),
0026 puppiPtCutsPhotons_(vd2vf(iConfig.getParameter<std::vector<double>>("puppiPtCutsPhotons"))),
0027 puppiUsingBareTracks_(iConfig.getParameter<bool>("puppiUsingBareTracks")) {
0028 debug_ = iConfig.getUntrackedParameter<int>("puppiDebug", debug_);
0029 if (puppiEtaCuts_.size() != puppiPtCuts_.size() || puppiPtCuts_.size() != puppiPtCutsPhotons_.size()) {
0030 throw cms::Exception("Configuration", "Bad PUPPI config");
0031 }
0032 for (unsigned int i = 0, n = puppiEtaCuts_.size(); i < n; ++i) {
0033 intPuppiEtaCuts_.push_back(std::round(puppiEtaCuts_[i] * CaloCluster::ETAPHI_SCALE));
0034 intPuppiPtCuts_.push_back(std::round(puppiPtCuts_[i] * CaloCluster::PT_SCALE));
0035 intPuppiPtCutsPhotons_.push_back(std::round(puppiPtCutsPhotons_[i] * CaloCluster::PT_SCALE));
0036 }
0037 }
0038
0039 PuppiAlgo::~PuppiAlgo() {}
0040
0041 const std::vector<std::string> &PuppiAlgo::puGlobalNames() const {
0042 static const std::vector<std::string> names_{"alphaCMed", "alphaCRms", "alphaFMed", "alphaFRms"};
0043 return names_;
0044 }
0045 void PuppiAlgo::doPUGlobals(const std::vector<Region> &rs, float z0, float npu, std::vector<float> &globals) const {
0046 globals.resize(4);
0047 computePuppiMedRMS(rs, globals[0], globals[1], globals[2], globals[3]);
0048 }
0049
0050 void PuppiAlgo::runNeutralsPU(Region &r, float z0, float npu, const std::vector<float> &globals) const {
0051 std::vector<float> alphaC, alphaF;
0052 computePuppiAlphas(r, alphaC, alphaF);
0053 computePuppiWeights(r, alphaC, alphaF, globals[0], globals[1], globals[2], globals[3]);
0054 fillPuppi(r);
0055 }
0056
0057 void PuppiAlgo::runNeutralsPU(Region &r, std::vector<float> &z0, float npu, const std::vector<float> &globals) const {
0058 float z0tmp = 0;
0059 runNeutralsPU(r, z0tmp, npu, globals);
0060 }
0061
0062 void PuppiAlgo::computePuppiAlphas(const Region &r, std::vector<float> &alphaC, std::vector<float> &alphaF) const {
0063 alphaC.resize(r.pf.size());
0064 alphaF.resize(r.pf.size());
0065 float puppiDr2 = std::pow(puppiDr_, 2), puppiDr2min = std::pow(puppiDrMin_, 2);
0066 for (unsigned int ip = 0, np = r.pf.size(); ip < np; ++ip) {
0067 const PFParticle &p = r.pf[ip];
0068 if (p.hwId <= 1)
0069 continue;
0070
0071 alphaC[ip] = 0;
0072 alphaF[ip] = 0;
0073 for (const PFParticle &p2 : r.pf) {
0074 float dr2 = ::deltaR2(p.floatEta(), p.floatPhi(), p2.floatEta(), p2.floatPhi());
0075 if (dr2 > 0 && dr2 < puppiDr2) {
0076 float w = std::pow(std::min(p2.floatPt(), puppiPtMax_), 2) / std::max<float>(puppiDr2min, dr2);
0077 alphaF[ip] += w;
0078 if (p2.chargedPV)
0079 alphaC[ip] += w;
0080 }
0081 }
0082 if (puppiUsingBareTracks_) {
0083 alphaC[ip] = 0;
0084 for (const PropagatedTrack &p2 : r.track) {
0085 if (!p2.fromPV)
0086 continue;
0087 if (!p2.quality(l1tpf_impl::InputTrack::PFLOOSE))
0088 continue;
0089 float dr2 = ::deltaR2(p.floatEta(), p.floatPhi(), p2.floatEta(), p2.floatPhi());
0090 if (dr2 > 0 && dr2 < puppiDr2) {
0091 alphaC[ip] += std::pow(std::min(p2.floatPt(), puppiPtMax_), 2) / std::max<float>(puppiDr2min, dr2);
0092 }
0093 }
0094 }
0095 }
0096 }
0097
0098 void PuppiAlgo::computePuppiWeights(Region &r,
0099 const std::vector<float> &alphaC,
0100 const std::vector<float> &alphaF,
0101 float alphaCMed,
0102 float alphaCRms,
0103 float alphaFMed,
0104 float alphaFRms) const {
0105 int16_t ietacut = std::round(etaCharged_ * CaloCluster::ETAPHI_SCALE);
0106 for (unsigned int ip = 0, np = r.pf.size(); ip < np; ++ip) {
0107 PFParticle &p = r.pf[ip];
0108
0109 if (p.hwId == l1t::PFCandidate::ChargedHadron || p.hwId == l1t::PFCandidate::Electron ||
0110 p.hwId == l1t::PFCandidate::Muon) {
0111 p.setPuppiW(p.chargedPV || p.hwId == l1t::PFCandidate::Muon ? 1.0 : 0);
0112 if (debug_)
0113 dbgPrintf(
0114 "PUPPI \t charged id %1d pt %7.2f eta %+5.2f phi %+5.2f alpha %+7.2f x2 %+7.2f --> puppi weight %.3f "
0115 "puppi pt %7.2f \n",
0116 p.hwId,
0117 p.floatPt(),
0118 p.floatEta(),
0119 p.floatPhi(),
0120 0.,
0121 0.,
0122 p.floatPuppiW(),
0123 p.floatPt() * p.floatPuppiW());
0124 continue;
0125 }
0126
0127 float alpha = -99, x2 = -99;
0128 bool central = std::abs(p.hwEta) < ietacut;
0129 if (r.relativeCoordinates)
0130 central =
0131 (std::abs(r.globalAbsEta(p.floatEta())) < etaCharged_);
0132 if (central) {
0133 if (alphaC[ip] > 0) {
0134 alpha = std::log(alphaC[ip]);
0135 x2 = (alpha - alphaCMed) * std::abs(alpha - alphaCMed) / std::pow(alphaCRms, 2);
0136 p.setPuppiW(ROOT::Math::chisquared_cdf(x2, 1));
0137 } else {
0138 p.setPuppiW(0);
0139 }
0140 } else {
0141 if (alphaF[ip] > 0) {
0142 alpha = std::log(alphaF[ip]);
0143 x2 = (alpha - alphaFMed) * std::abs(alpha - alphaFMed) / std::pow(alphaFRms, 2);
0144 p.setPuppiW(ROOT::Math::chisquared_cdf(x2, 1));
0145 } else {
0146 p.setPuppiW(0);
0147 }
0148 }
0149 if (debug_)
0150 dbgPrintf(
0151 "PUPPI \t neutral id %1d pt %7.2f eta %+5.2f phi %+5.2f alpha %+7.2f x2 %+7.2f --> puppi weight %.3f "
0152 "puppi pt %7.2f \n",
0153 p.hwId,
0154 p.floatPt(),
0155 p.floatEta(),
0156 p.floatPhi(),
0157 alpha,
0158 x2,
0159 p.floatPuppiW(),
0160 p.floatPt() * p.floatPuppiW());
0161 }
0162 }
0163
0164 void PuppiAlgo::computePuppiMedRMS(
0165 const std::vector<Region> &rs, float &alphaCMed, float &alphaCRms, float &alphaFMed, float &alphaFRms) const {
0166 std::vector<float> alphaFs;
0167 std::vector<float> alphaCs;
0168 int16_t ietacut = std::round(etaCharged_ * CaloCluster::ETAPHI_SCALE);
0169 float puppiDr2 = std::pow(puppiDr_, 2), puppiDr2min = std::pow(puppiDrMin_, 2);
0170 for (const Region &r : rs) {
0171 for (const PFParticle &p : r.pf) {
0172 bool central = std::abs(p.hwEta) < ietacut;
0173 if (r.relativeCoordinates)
0174 central = (r.globalAbsEta(p.floatEta()) < etaCharged_);
0175 if (central) {
0176 if (p.hwId > 1 || p.chargedPV)
0177 continue;
0178 }
0179 float alphaC = 0, alphaF = 0;
0180 for (const PFParticle &p2 : r.pf) {
0181 float dr2 = ::deltaR2(p.floatEta(), p.floatPhi(), p2.floatEta(), p2.floatPhi());
0182 if (dr2 > 0 && dr2 < puppiDr2) {
0183 float w = std::pow(std::min(p2.floatPt(), puppiPtMax_), 2) / std::max<float>(puppiDr2min, dr2);
0184 alphaF += w;
0185 if (p2.chargedPV)
0186 alphaC += w;
0187 }
0188 }
0189 if (puppiUsingBareTracks_) {
0190 alphaC = 0;
0191 for (const PropagatedTrack &p2 : r.track) {
0192 if (!p2.fromPV)
0193 continue;
0194 float dr2 = ::deltaR2(p.floatEta(), p.floatPhi(), p2.floatEta(), p2.floatPhi());
0195 if (dr2 > 0 && dr2 < puppiDr2) {
0196 alphaC += std::pow(std::min(p2.floatPt(), puppiPtMax_), 2) / std::max<float>(puppiDr2min, dr2);
0197 }
0198 }
0199 }
0200 if (central) {
0201 if (alphaC > 0)
0202 alphaCs.push_back(std::log(alphaC));
0203 } else {
0204 if (alphaF > 0)
0205 alphaFs.push_back(std::log(alphaF));
0206 }
0207 }
0208 }
0209 std::sort(alphaCs.begin(), alphaCs.end());
0210 std::sort(alphaFs.begin(), alphaFs.end());
0211
0212 if (alphaCs.size() > 1) {
0213 alphaCMed = alphaCs[alphaCs.size() / 2 + 1];
0214 double sum = 0.0;
0215 for (float alpha : alphaCs)
0216 sum += std::pow(alpha - alphaCMed, 2);
0217 alphaCRms = std::sqrt(float(sum) / alphaCs.size());
0218 } else {
0219 alphaCMed = 8.;
0220 alphaCRms = 8.;
0221 }
0222
0223 if (alphaFs.size() > 1) {
0224 alphaFMed = alphaFs[alphaFs.size() / 2 + 1];
0225 double sum = 0.0;
0226 for (float alpha : alphaFs)
0227 sum += std::pow(alpha - alphaFMed, 2);
0228 alphaFRms = std::sqrt(float(sum) / alphaFs.size());
0229 } else {
0230 alphaFMed = 6.;
0231 alphaFRms = 6.;
0232 }
0233 if (debug_)
0234 dbgPrintf("PUPPI \t alphaC = %+6.2f +- %6.2f (%4lu), alphaF = %+6.2f +- %6.2f (%4lu)\n",
0235 alphaCMed,
0236 alphaCRms,
0237 alphaCs.size(),
0238 alphaFMed,
0239 alphaFRms,
0240 alphaFs.size());
0241 }
0242
0243 void PuppiAlgo::fillPuppi(Region &r) const {
0244 uint16_t PUPPIW_0p01 = std::round(0.01 * PFParticle::PUPPI_SCALE);
0245 r.puppi.clear();
0246 for (PFParticle &p : r.pf) {
0247 if (p.hwId == l1t::PFCandidate::ChargedHadron || p.hwId == l1t::PFCandidate::Electron ||
0248 p.hwId == l1t::PFCandidate::Muon) {
0249 if (p.hwPuppiWeight > 0) {
0250 r.puppi.push_back(p);
0251 }
0252 } else {
0253 if (p.hwPuppiWeight > PUPPIW_0p01) {
0254
0255
0256 int16_t hwPt = (float(p.hwPt) * float(p.hwPuppiWeight) / float(PFParticle::PUPPI_SCALE));
0257 int16_t hwPtCut = 0, hwAbsEta = r.relativeCoordinates
0258 ? round(r.globalAbsEta(p.floatEta()) * CaloCluster::ETAPHI_SCALE)
0259 : std::abs(p.hwEta);
0260 for (unsigned int ietaBin = 0, nBins = intPuppiEtaCuts_.size(); ietaBin < nBins; ++ietaBin) {
0261 if (hwAbsEta < intPuppiEtaCuts_[ietaBin]) {
0262 hwPtCut = (p.hwId == l1t::PFCandidate::Photon ? intPuppiPtCutsPhotons_[ietaBin] : intPuppiPtCuts_[ietaBin]);
0263 break;
0264 }
0265 }
0266 if (hwPt > hwPtCut) {
0267 r.puppi.push_back(p);
0268 r.puppi.back().hwPt = hwPt;
0269 }
0270 }
0271 }
0272 }
0273 }