File indexing completed on 2024-09-07 04:35:57
0001 #ifndef RecoCandidate_IsoDeposit_H
0002 #define RecoCandidate_IsoDeposit_H
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019 #include "DataFormats/RecoCandidate/interface/IsoDepositDirection.h"
0020 #include "DataFormats/Math/interface/Vector3D.h"
0021 #include "FWCore/Utilities/interface/Exception.h"
0022 #include "FWCore/Utilities/interface/thread_safety_macros.h"
0023 #include <map>
0024 #include <cmath>
0025 #include <string>
0026 #include <vector>
0027 #include <typeinfo>
0028 #include <atomic>
0029
0030 namespace reco {
0031 namespace isodeposit {
0032 struct AbsVeto {
0033 virtual ~AbsVeto() {}
0034
0035 virtual bool veto(double eta, double phi, float value) const = 0;
0036
0037
0038 virtual void centerOn(double eta, double phi) {
0039 throw cms::Exception("Not Implemented") << "This AbsVeto implementation (" << typeid(this).name()
0040 << ") does not support the centerOn(eta,phi) method";
0041 }
0042 };
0043 typedef std::vector<AbsVeto*> AbsVetos;
0044 }
0045 }
0046
0047 namespace reco {
0048
0049 class IsoDeposit {
0050 public:
0051 typedef isodeposit::Direction Direction;
0052 typedef isodeposit::AbsVeto AbsVeto;
0053 typedef isodeposit::AbsVetos AbsVetos;
0054 typedef Direction::Distance Distance;
0055 typedef std::multimap<Distance, float> DepositsMultimap;
0056 typedef DepositsMultimap::const_iterator DepIterator;
0057
0058
0059 struct Veto {
0060 Direction vetoDir;
0061 float dR;
0062 Veto() {}
0063 Veto(Direction dir, double d) : vetoDir(dir), dR(d) {}
0064 };
0065 typedef std::vector<Veto> Vetos;
0066
0067
0068 IsoDeposit(double eta = 0, double phi = 0);
0069 IsoDeposit(const Direction& candDirection);
0070
0071
0072 virtual ~IsoDeposit() {}
0073
0074
0075 const Direction& direction() const { return theDirection; }
0076 double eta() const { return theDirection.eta(); }
0077 double phi() const { return theDirection.phi(); }
0078
0079
0080 const Veto& veto() const { return theVeto; }
0081
0082 void setVeto(const Veto& aVeto) { theVeto = aVeto; }
0083
0084
0085 void addDeposit(double dr, double deposit);
0086 void addDeposit(const Direction& depDir, double deposit);
0087
0088
0089 double depositWithin(double coneSize,
0090 const Vetos& vetos = Vetos(),
0091 bool skipDepositVeto = false
0092 ) const;
0093
0094
0095 double depositWithin(Direction dir,
0096 double coneSize,
0097 const Vetos& vetos = Vetos(),
0098 bool skipDepositVeto = false
0099 ) const;
0100
0101
0102 std::pair<double, int> depositAndCountWithin(double coneSize,
0103 const Vetos& vetos = Vetos(),
0104 double threshold = -1e+36,
0105 bool skipDepositVeto = false
0106 ) const;
0107
0108
0109 std::pair<double, int> depositAndCountWithin(Direction dir,
0110 double coneSize,
0111 const Vetos& vetos = Vetos(),
0112 double threshold = -1e+36,
0113 bool skipDepositVeto = false
0114 ) const;
0115
0116
0117 double depositWithin(double coneSize,
0118 const AbsVetos& vetos,
0119 bool skipDepositVeto = false
0120 ) const;
0121
0122
0123 std::pair<double, int> depositAndCountWithin(double coneSize,
0124 const AbsVetos& vetos,
0125 bool skipDepositVeto = false
0126 ) const;
0127
0128
0129 double candEnergy() const { return theCandTag; }
0130
0131
0132 void addCandEnergy(double et) { theCandTag += et; }
0133
0134 std::string print() const;
0135
0136 class const_iterator {
0137 public:
0138 const const_iterator& operator++() {
0139 ++it_;
0140 cacheReady_ = false;
0141 return *this;
0142 }
0143 const const_iterator* operator->() const { return this; }
0144 float dR() const { return it_->first.deltaR; }
0145 float eta() const {
0146 if (!cacheReady_)
0147 doDir();
0148 return cache_.eta();
0149 }
0150 float phi() const {
0151 if (!cacheReady_)
0152 doDir();
0153 return cache_.phi();
0154 }
0155 float value() const { return it_->second; }
0156 bool operator!=(const const_iterator& it2) { return it2.it_ != it_; }
0157 friend class IsoDeposit;
0158
0159 private:
0160 typedef Direction::Distance Distance;
0161 void doDir() const {
0162 if (!cacheReady_)
0163 cache_ = parent_->direction() + it_->first;
0164 cacheReady_ = true;
0165 }
0166 const_iterator(const IsoDeposit* parent, std::multimap<Distance, float>::const_iterator it)
0167 : parent_(parent), it_(it), cache_(), cacheReady_(false) {}
0168 const reco::IsoDeposit* parent_;
0169 std::multimap<Distance, float>::const_iterator it_;
0170 CMS_THREAD_SAFE mutable Direction cache_;
0171 mutable std::atomic<bool> cacheReady_;
0172 };
0173 const_iterator begin() const { return const_iterator(this, theDeposits.begin()); }
0174 const_iterator end() const { return const_iterator(this, theDeposits.end()); }
0175
0176 class SumAlgo {
0177 public:
0178 SumAlgo() : sum_(0) {}
0179 void operator+=(DepIterator deposit) { sum_ += deposit->second; }
0180 void operator+=(double deposit) { sum_ += deposit; }
0181 double result() const { return sum_; }
0182
0183 private:
0184 double sum_;
0185 };
0186 class CountAlgo {
0187 public:
0188 CountAlgo() : count_(0) {}
0189 void operator+=(DepIterator deposit) { count_++; }
0190 void operator+=(double deposit) { count_++; }
0191 double result() const { return count_; }
0192
0193 private:
0194 size_t count_;
0195 };
0196 class Sum2Algo {
0197 public:
0198 Sum2Algo() : sum2_(0) {}
0199 void operator+=(DepIterator deposit) { sum2_ += deposit->second * deposit->second; }
0200 void operator+=(double deposit) { sum2_ += deposit * deposit; }
0201 double result() const { return sum2_; }
0202
0203 private:
0204 double sum2_;
0205 };
0206 class MaxAlgo {
0207 public:
0208 MaxAlgo() : max_(0) {}
0209 void operator+=(DepIterator deposit) {
0210 if (deposit->second > max_)
0211 max_ = deposit->second;
0212 }
0213 void operator+=(double deposit) {
0214 if (deposit > max_)
0215 max_ = deposit;
0216 }
0217 double result() const { return max_; }
0218
0219 private:
0220 double max_;
0221 };
0222
0223 class MeanDRAlgo {
0224 public:
0225 MeanDRAlgo() : sum_(0.), count_(0.) {}
0226 void operator+=(DepIterator deposit) {
0227 sum_ += deposit->first.deltaR;
0228 count_ += 1.0;
0229 }
0230 double result() const { return sum_ / std::max(1., count_); }
0231
0232 private:
0233 double sum_;
0234 double count_;
0235 };
0236
0237 class SumDRAlgo {
0238 public:
0239 SumDRAlgo() : sum_(0.) {}
0240 void operator+=(DepIterator deposit) { sum_ += deposit->first.deltaR; }
0241 double result() const { return sum_; }
0242
0243 private:
0244 double sum_;
0245 };
0246
0247
0248 template <typename Algo>
0249 double algoWithin(double coneSize,
0250 const AbsVetos& vetos = AbsVetos(),
0251 bool skipDepositVeto = false
0252 ) const;
0253
0254 template <typename Algo>
0255 double algoWithin(const Direction&,
0256 double coneSize,
0257 const AbsVetos& vetos = AbsVetos(),
0258 bool skipDepositVeto = false
0259 ) const;
0260
0261 double countWithin(double coneSize,
0262 const AbsVetos& vetos = AbsVetos(),
0263 bool skipDepositVeto = false
0264 ) const;
0265
0266 double sumWithin(double coneSize,
0267 const AbsVetos& vetos = AbsVetos(),
0268 bool skipDepositVeto = false
0269 ) const;
0270
0271 double sumWithin(const Direction& dir,
0272 double coneSize,
0273 const AbsVetos& vetos = AbsVetos(),
0274 bool skipDepositVeto = false
0275 ) const;
0276 double sum2Within(double coneSize,
0277 const AbsVetos& vetos = AbsVetos(),
0278 bool skipDepositVeto = false
0279 ) const;
0280
0281 double maxWithin(double coneSize,
0282 const AbsVetos& vetos = AbsVetos(),
0283 bool skipDepositVeto = false
0284 ) const;
0285
0286 double nearestDR(double coneSize,
0287 const AbsVetos& vetos = AbsVetos(),
0288 bool skipDepositVeto = false
0289 ) const;
0290
0291 private:
0292
0293 Direction theDirection;
0294
0295
0296 Veto theVeto;
0297
0298
0299 float theCandTag;
0300
0301
0302
0303 DepositsMultimap theDeposits;
0304 };
0305
0306 }
0307
0308 template <typename Algo>
0309 double reco::IsoDeposit::algoWithin(double coneSize, const AbsVetos& vetos, bool skipDepositVeto) const {
0310 using namespace reco::isodeposit;
0311 Algo algo;
0312 typedef AbsVetos::const_iterator IV;
0313 IV ivEnd = vetos.end();
0314
0315 Distance maxDistance = {float(coneSize), 999.f};
0316 typedef DepositsMultimap::const_iterator IM;
0317 IM imLoc = theDeposits.upper_bound(maxDistance);
0318 for (IM im = theDeposits.begin(); im != imLoc; ++im) {
0319 bool vetoed = false;
0320 Direction dirDep = theDirection + im->first;
0321 for (IV iv = vetos.begin(); iv < ivEnd; ++iv) {
0322 if ((*iv)->veto(dirDep.eta(), dirDep.phi(), im->second)) {
0323 vetoed = true;
0324 break;
0325 }
0326 }
0327 if (!vetoed) {
0328 if (skipDepositVeto || (dirDep.deltaR(theVeto.vetoDir) > theVeto.dR)) {
0329 algo += im;
0330 }
0331 }
0332 }
0333 return algo.result();
0334 }
0335
0336 template <typename Algo>
0337 double reco::IsoDeposit::algoWithin(const Direction& dir,
0338 double coneSize,
0339 const AbsVetos& vetos,
0340 bool skipDepositVeto) const {
0341 using namespace reco::isodeposit;
0342 Algo algo;
0343 typedef AbsVetos::const_iterator IV;
0344 IV ivEnd = vetos.end();
0345 typedef DepositsMultimap::const_iterator IM;
0346 IM imLoc = theDeposits.end();
0347 for (IM im = theDeposits.begin(); im != imLoc; ++im) {
0348 bool vetoed = false;
0349 Direction dirDep = theDirection + im->first;
0350 Distance newDist = dirDep - dir;
0351 if (newDist.deltaR > coneSize)
0352 continue;
0353 for (IV iv = vetos.begin(); iv < ivEnd; ++iv) {
0354 if ((*iv)->veto(dirDep.eta(), dirDep.phi(), im->second)) {
0355 vetoed = true;
0356 break;
0357 }
0358 }
0359 if (!vetoed) {
0360 if (skipDepositVeto || (dirDep.deltaR(theVeto.vetoDir) > theVeto.dR)) {
0361 algo += im;
0362 }
0363 }
0364 }
0365 return algo.result();
0366 }
0367
0368 #endif