File indexing completed on 2022-02-04 00:11:04
0001 #ifndef DataFormats_L1Trigger_HGCalClusterT_h
0002 #define DataFormats_L1Trigger_HGCalClusterT_h
0003
0004
0005 #include "DataFormats/Common/interface/Ptr.h"
0006 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0007 #include "DataFormats/L1Trigger/interface/L1Candidate.h"
0008 #include "DataFormats/L1THGCal/interface/HGCalTriggerCell.h"
0009 #include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h"
0010 #include "DataFormats/ForwardDetId/interface/HGCalTriggerDetId.h"
0011 #include "DataFormats/ForwardDetId/interface/HFNoseTriggerDetId.h"
0012
0013
0014 #include "Math/Vector3D.h"
0015
0016 #include <unordered_map>
0017
0018 namespace l1t {
0019 template <class C>
0020 class HGCalClusterT : public L1Candidate {
0021 public:
0022 typedef typename std::unordered_map<uint32_t, edm::Ptr<C>>::const_iterator const_iterator;
0023
0024 public:
0025 HGCalClusterT() {}
0026 HGCalClusterT(const LorentzVector p4, int pt = 0, int eta = 0, int phi = 0)
0027 : L1Candidate(p4, pt, eta, phi),
0028 valid_(true),
0029 detId_(0),
0030 centre_(0, 0, 0),
0031 centreProj_(0., 0., 0.),
0032 mipPt_(0),
0033 seedMipPt_(0),
0034 sumPt_(0) {}
0035
0036 HGCalClusterT(const edm::Ptr<C>& c, float fraction = 1.)
0037 : valid_(true),
0038 detId_(c->detId()),
0039 centre_(0., 0., 0.),
0040 centreProj_(0., 0., 0.),
0041 mipPt_(0.),
0042 seedMipPt_(0.),
0043 sumPt_(0.) {
0044 addConstituent(c, true, fraction);
0045 }
0046
0047 ~HGCalClusterT() override{};
0048
0049 const std::unordered_map<uint32_t, edm::Ptr<C>>& constituents() const { return constituents_; }
0050 const_iterator constituents_begin() const { return constituents_.begin(); }
0051 const_iterator constituents_end() const { return constituents_.end(); }
0052 unsigned size() const { return constituents_.size(); }
0053
0054 void addConstituent(const edm::Ptr<C>& c, bool updateCentre = true, float fraction = 1.) {
0055 double cMipt = c->mipPt() * fraction;
0056
0057 if (constituents_.empty()) {
0058 detId_ = DetId(c->detId());
0059 seedMipPt_ = cMipt;
0060
0061
0062 if (!updateCentre) {
0063 centre_ = c->position();
0064 }
0065 }
0066 updateP4AndPosition(c, updateCentre, fraction);
0067
0068 constituents_.emplace(c->detId(), c);
0069 constituentsFraction_.emplace(c->detId(), fraction);
0070 }
0071
0072 void removeConstituent(const edm::Ptr<C>& c, bool updateCentre = true) {
0073
0074 double fraction = 0;
0075 const auto& constituent_itr = constituents_.find(c->detId());
0076 const auto& fraction_itr = constituentsFraction_.find(c->detId());
0077 if (constituent_itr != constituents_.end()) {
0078
0079 fraction = fraction_itr->second;
0080 constituents_.erase(constituent_itr);
0081 constituentsFraction_.erase(fraction_itr);
0082
0083 updateP4AndPosition(c, updateCentre, -fraction);
0084 }
0085 }
0086
0087 bool valid() const { return valid_; }
0088 void setValid(bool valid) { valid_ = valid; }
0089
0090 double mipPt() const { return mipPt_; }
0091 double seedMipPt() const { return seedMipPt_; }
0092 uint32_t detId() const { return detId_.rawId(); }
0093 void setDetId(uint32_t id) { detId_ = id; }
0094 void setPt(double pt) { setP4(math::PtEtaPhiMLorentzVector(pt, eta(), phi(), mass())); }
0095 double sumPt() const { return sumPt_; }
0096
0097 double distance(const l1t::HGCalTriggerCell& tc) const { return (tc.position() - centre_).mag(); }
0098
0099 const GlobalPoint& position() const { return centre_; }
0100 const GlobalPoint& centre() const { return centre_; }
0101 const GlobalPoint& centreProj() const { return centreProj_; }
0102
0103 double hOverE() const {
0104 double pt_em = 0.;
0105 double pt_had = 0.;
0106 double hOe = 0.;
0107
0108 for (const auto& id_constituent : constituents()) {
0109 DetId id(id_constituent.first);
0110 auto id_fraction = constituentsFraction_.find(id_constituent.first);
0111 double fraction = (id_fraction != constituentsFraction_.end() ? id_fraction->second : 1.);
0112 if ((id.det() == DetId::HGCalEE) ||
0113 (id.det() == DetId::HGCalTrigger &&
0114 HGCalTriggerDetId(id).subdet() == HGCalTriggerSubdetector::HGCalEETrigger) ||
0115 (id.det() == DetId::Forward && id.subdetId() == ForwardSubdetector::HFNose && HFNoseDetId(id).isEE()) ||
0116 (id.det() == DetId::HGCalTrigger &&
0117 HGCalTriggerDetId(id).subdet() == HGCalTriggerSubdetector::HFNoseTrigger &&
0118 HFNoseTriggerDetId(id).isEE())) {
0119 pt_em += id_constituent.second->pt() * fraction;
0120 } else if ((id.det() == DetId::HGCalHSi) || (id.det() == DetId::HGCalHSc) ||
0121 (id.det() == DetId::HGCalTrigger &&
0122 HGCalTriggerDetId(id).subdet() == HGCalTriggerSubdetector::HGCalHSiTrigger) ||
0123 (id.det() == DetId::Forward && id.subdetId() == ForwardSubdetector::HFNose &&
0124 HFNoseDetId(id).isHE()) ||
0125 (id.det() == DetId::HGCalTrigger &&
0126 HGCalTriggerDetId(id).subdet() == HGCalTriggerSubdetector::HFNoseTrigger &&
0127 HFNoseTriggerDetId(id).isHSilicon())) {
0128 pt_had += id_constituent.second->pt() * fraction;
0129 }
0130 }
0131 if (pt_em > 0)
0132 hOe = pt_had / pt_em;
0133 else
0134 hOe = -1.;
0135 return hOe;
0136 }
0137
0138 uint32_t subdetId() const { return detId_.subdetId(); }
0139
0140
0141
0142 int showerLength() const { return showerLength_; }
0143 int coreShowerLength() const { return coreShowerLength_; }
0144 int firstLayer() const { return firstLayer_; }
0145 int maxLayer() const { return maxLayer_; }
0146 float eMax() const { return eMax_; }
0147 float sigmaEtaEtaMax() const { return sigmaEtaEtaMax_; }
0148 float sigmaPhiPhiMax() const { return sigmaPhiPhiMax_; }
0149 float sigmaEtaEtaTot() const { return sigmaEtaEtaTot_; }
0150 float sigmaPhiPhiTot() const { return sigmaPhiPhiTot_; }
0151 float sigmaZZ() const { return sigmaZZ_; }
0152 float sigmaRRTot() const { return sigmaRRTot_; }
0153 float sigmaRRMax() const { return sigmaRRMax_; }
0154 float sigmaRRMean() const { return sigmaRRMean_; }
0155 float zBarycenter() const { return zBarycenter_; }
0156 float layer10percent() const { return layer10percent_; }
0157 float layer50percent() const { return layer50percent_; }
0158 float layer90percent() const { return layer90percent_; }
0159 float triggerCells67percent() const { return triggerCells67percent_; }
0160 float triggerCells90percent() const { return triggerCells90percent_; }
0161
0162 void showerLength(int showerLength) { showerLength_ = showerLength; }
0163 void coreShowerLength(int coreShowerLength) { coreShowerLength_ = coreShowerLength; }
0164 void firstLayer(int firstLayer) { firstLayer_ = firstLayer; }
0165 void maxLayer(int maxLayer) { maxLayer_ = maxLayer; }
0166 void eMax(float eMax) { eMax_ = eMax; }
0167 void sigmaEtaEtaMax(float sigmaEtaEtaMax) { sigmaEtaEtaMax_ = sigmaEtaEtaMax; }
0168 void sigmaEtaEtaTot(float sigmaEtaEtaTot) { sigmaEtaEtaTot_ = sigmaEtaEtaTot; }
0169 void sigmaPhiPhiMax(float sigmaPhiPhiMax) { sigmaPhiPhiMax_ = sigmaPhiPhiMax; }
0170 void sigmaPhiPhiTot(float sigmaPhiPhiTot) { sigmaPhiPhiTot_ = sigmaPhiPhiTot; }
0171 void sigmaRRMax(float sigmaRRMax) { sigmaRRMax_ = sigmaRRMax; }
0172 void sigmaRRTot(float sigmaRRTot) { sigmaRRTot_ = sigmaRRTot; }
0173 void sigmaRRMean(float sigmaRRMean) { sigmaRRMean_ = sigmaRRMean; }
0174 void sigmaZZ(float sigmaZZ) { sigmaZZ_ = sigmaZZ; }
0175 void zBarycenter(float zBarycenter) { zBarycenter_ = zBarycenter; }
0176 void layer10percent(float layer10percent) { layer10percent_ = layer10percent; }
0177 void layer50percent(float layer50percent) { layer50percent_ = layer50percent; }
0178 void layer90percent(float layer90percent) { layer90percent_ = layer90percent; }
0179 void triggerCells67percent(float triggerCells67percent) { triggerCells67percent_ = triggerCells67percent; }
0180 void triggerCells90percent(float triggerCells90percent) { triggerCells90percent_ = triggerCells90percent; }
0181
0182
0183 bool operator<(const HGCalClusterT<C>& cl) const { return mipPt() < cl.mipPt(); }
0184 bool operator>(const HGCalClusterT<C>& cl) const { return cl < *this; }
0185 bool operator<=(const HGCalClusterT<C>& cl) const { return !(cl > *this); }
0186 bool operator>=(const HGCalClusterT<C>& cl) const { return !(cl < *this); }
0187
0188 private:
0189 bool valid_ = false;
0190 DetId detId_;
0191
0192 std::unordered_map<uint32_t, edm::Ptr<C>> constituents_;
0193 std::unordered_map<uint32_t, double> constituentsFraction_;
0194
0195 GlobalPoint centre_;
0196 GlobalPoint centreProj_;
0197
0198 double mipPt_ = 0.;
0199 double seedMipPt_ = 0.;
0200 double sumPt_ = 0.;
0201
0202
0203
0204 int showerLength_ = 0;
0205 int coreShowerLength_ = 0;
0206 int firstLayer_ = 0;
0207 int maxLayer_ = 0;
0208 float eMax_ = 0.;
0209 float sigmaEtaEtaMax_ = 0.;
0210 float sigmaPhiPhiMax_ = 0.;
0211 float sigmaRRMax_ = 0.;
0212 float sigmaEtaEtaTot_ = 0.;
0213 float sigmaPhiPhiTot_ = 0.;
0214 float sigmaRRTot_ = 0.;
0215 float sigmaRRMean_ = 0.;
0216 float sigmaZZ_ = 0.;
0217 float zBarycenter_ = 0.;
0218 float layer10percent_ = 0.;
0219 float layer50percent_ = 0.;
0220 float layer90percent_ = 0.;
0221 float triggerCells67percent_ = 0.;
0222 float triggerCells90percent_ = 0.;
0223
0224 void updateP4AndPosition(const edm::Ptr<C>& c, bool updateCentre = true, float fraction = 1.) {
0225 double cMipt = c->mipPt() * fraction;
0226 double cPt = c->pt() * fraction;
0227
0228 if (updateCentre) {
0229 Basic3DVector<float> constituentCentre(c->position());
0230 Basic3DVector<float> clusterCentre(centre_);
0231
0232 clusterCentre = clusterCentre * mipPt_ + constituentCentre * cMipt;
0233 if ((mipPt_ + cMipt) > 0) {
0234 clusterCentre /= (mipPt_ + cMipt);
0235 }
0236 centre_ = GlobalPoint(clusterCentre);
0237
0238 if (clusterCentre.z() != 0) {
0239 centreProj_ = GlobalPoint(clusterCentre / std::abs(clusterCentre.z()));
0240 }
0241 }
0242
0243
0244 mipPt_ += cMipt;
0245 sumPt_ += cPt;
0246 int updatedPt = hwPt() + (int)(c->hwPt() * fraction);
0247 setHwPt(updatedPt);
0248
0249 math::PtEtaPhiMLorentzVector updatedP4(p4());
0250 updatedP4 += (c->p4() * fraction);
0251 setP4(updatedP4);
0252 }
0253 };
0254
0255 }
0256
0257 #endif