Back to home page

Project CMSSW displayed by LXR
 
 

    

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 }  // namespace
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     // neutral
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     // charged
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     // neutral
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_);  // FIXME could make a better integer implementation
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_);  // FIXME could make a better integer implementation
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) {  // charged
0249       if (p.hwPuppiWeight > 0) {
0250         r.puppi.push_back(p);
0251       }
0252     } else {  // neutral
0253       if (p.hwPuppiWeight > PUPPIW_0p01) {
0254         // FIXME would work better with PUPPI_SCALE being a power of two, to do the shift
0255         // FIXME done with floats
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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!