L1THGCal/interface/HGCalClusterT.h

0001 #ifndef DataFormats_L1Trigger_HGCalClusterT_h
0002 #define DataFormats_L1Trigger_HGCalClusterT_h
0003
0004 /* CMSSW */
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 /* ROOT */
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         /* if the centre will not be dynamically calculated
0061              the seed centre is considere as cluster centre */
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       /* remove the pointer to c from the std::vector */
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         // remove constituent and get its fraction in the cluster
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     /* distance in 'cm' */
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     //shower shape
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 varRR() const { return varRR_; }
0154     float varZZ() const { return varZZ_; }
0155     float varEtaEta() const { return varEtaEta_; }
0156     float varPhiPhi() const { return varPhiPhi_; }
0157     float sigmaRRMax() const { return sigmaRRMax_; }
0158     float sigmaRRMean() const { return sigmaRRMean_; }
0159     float zBarycenter() const { return zBarycenter_; }
0160     float layer10percent() const { return layer10percent_; }
0161     float layer50percent() const { return layer50percent_; }
0162     float layer90percent() const { return layer90percent_; }
0163     float triggerCells67percent() const { return triggerCells67percent_; }
0164     float triggerCells90percent() const { return triggerCells90percent_; }
0165     float first1layers() const { return first1layers_; }
0166     float first3layers() const { return first3layers_; }
0167     float first5layers() const { return first5layers_; }
0168     float firstHcal1layers() const { return firstHcal1layers_; }
0169     float firstHcal3layers() const { return firstHcal3layers_; }
0170     float firstHcal5layers() const { return firstHcal5layers_; }
0171     float last1layers() const { return last1layers_; }
0172     float last3layers() const { return last3layers_; }
0173     float last5layers() const { return last5layers_; }
0174     float emax1layers() const { return emax1layers_; }
0175     float emax3layers() const { return emax3layers_; }
0176     float emax5layers() const { return emax5layers_; }
0177     float eot() const { return eot_; }
0178     int ebm0() const { return ebm0_; }
0179     int ebm1() const { return ebm1_; }
0180     int hbm() const { return hbm_; }
0181
0182     void setShowerLength(int showerLength) { showerLength_ = showerLength; }
0183     void setCoreShowerLength(int coreShowerLength) { coreShowerLength_ = coreShowerLength; }
0184     void setFirstLayer(int firstLayer) { firstLayer_ = firstLayer; }
0185     void setMaxLayer(int maxLayer) { maxLayer_ = maxLayer; }
0186     void setEMax(float eMax) { eMax_ = eMax; }
0187     void setSigmaEtaEtaMax(float sigmaEtaEtaMax) { sigmaEtaEtaMax_ = sigmaEtaEtaMax; }
0188     void setSigmaEtaEtaTot(float sigmaEtaEtaTot) { sigmaEtaEtaTot_ = sigmaEtaEtaTot; }
0189     void setSigmaPhiPhiMax(float sigmaPhiPhiMax) { sigmaPhiPhiMax_ = sigmaPhiPhiMax; }
0190     void setSigmaPhiPhiTot(float sigmaPhiPhiTot) { sigmaPhiPhiTot_ = sigmaPhiPhiTot; }
0191     void setSigmaRRMax(float sigmaRRMax) { sigmaRRMax_ = sigmaRRMax; }
0192     void setSigmaRRTot(float sigmaRRTot) { sigmaRRTot_ = sigmaRRTot; }
0193     void setVarRR(float varRR) { varRR_ = varRR; }
0194     void setVarZZ(float varZZ) { varZZ_ = varZZ; }
0195     void setVarEtaEta(float varEtaEta) { varEtaEta_ = varEtaEta; }
0196     void setVarPhiPhi(float varPhiPhi) { varPhiPhi_ = varPhiPhi; }
0197     void setSigmaRRMean(float sigmaRRMean) { sigmaRRMean_ = sigmaRRMean; }
0198     void setSigmaZZ(float sigmaZZ) { sigmaZZ_ = sigmaZZ; }
0199     void setZBarycenter(float zBarycenter) { zBarycenter_ = zBarycenter; }
0200     void setLayer10percent(float layer10percent) { layer10percent_ = layer10percent; }
0201     void setLayer50percent(float layer50percent) { layer50percent_ = layer50percent; }
0202     void setLayer90percent(float layer90percent) { layer90percent_ = layer90percent; }
0203     void setTriggerCells67percent(float triggerCells67percent) { triggerCells67percent_ = triggerCells67percent; }
0204     void setTriggerCells90percent(float triggerCells90percent) { triggerCells90percent_ = triggerCells90percent; }
0205     void setFirst1layers(float first1layers) { first1layers_ = first1layers; }
0206     void setFirst3layers(float first3layers) { first3layers_ = first3layers; }
0207     void setFirst5layers(float first5layers) { first5layers_ = first5layers; }
0208     void setFirstHcal1layers(float firstHcal1layers) { firstHcal1layers_ = firstHcal1layers; }
0209     void setFirstHcal3layers(float firstHcal3layers) { firstHcal3layers_ = firstHcal3layers; }
0210     void setFirstHcal5layers(float firstHcal5layers) { firstHcal5layers_ = firstHcal5layers; }
0211     void setLast1layers(float last1layers) { last1layers_ = last1layers; }
0212     void setLast3layers(float last3layers) { last3layers_ = last3layers; }
0213     void setLast5layers(float last5layers) { last5layers_ = last5layers; }
0214     void setEmax1layers(float emax1layers) { emax1layers_ = emax1layers; }
0215     void setEmax3layers(float emax3layers) { emax3layers_ = emax3layers; }
0216     void setEmax5layers(float emax5layers) { emax5layers_ = emax5layers; }
0217     void setEot(float eot) { eot_ = eot; }
0218     void setEbm0(int ebm0) { ebm0_ = ebm0; }
0219     void setEbm1(int ebm1) { ebm1_ = ebm1; }
0220     void setHbm(int hbm) { hbm_ = hbm; }
0221
0222     /* operators */
0223     bool operator<(const HGCalClusterT<C>& cl) const { return mipPt() < cl.mipPt(); }
0224     bool operator>(const HGCalClusterT<C>& cl) const { return cl < *this; }
0225     bool operator<=(const HGCalClusterT<C>& cl) const { return !(cl > *this); }
0226     bool operator>=(const HGCalClusterT<C>& cl) const { return !(cl < *this); }
0227
0228   private:
0229     bool valid_ = false;
0230     DetId detId_;
0231
0232     std::unordered_map<uint32_t, edm::Ptr<C>> constituents_;
0233     std::unordered_map<uint32_t, double> constituentsFraction_;
0234
0235     GlobalPoint centre_;
0236     GlobalPoint centreProj_;  // centre projected onto the first HGCal layer
0237
0238     double mipPt_ = 0.;
0239     double seedMipPt_ = 0.;
0240     double sumPt_ = 0.;
0241
0242     //shower shape
0243
0244     int showerLength_ = 0;
0245     int coreShowerLength_ = 0;
0246     int firstLayer_ = 0;
0247     int maxLayer_ = 0;
0248     float eMax_ = 0.;
0249     float sigmaEtaEtaMax_ = 0.;
0250     float sigmaPhiPhiMax_ = 0.;
0251     float sigmaRRMax_ = 0.;
0252     float sigmaEtaEtaTot_ = 0.;
0253     float sigmaPhiPhiTot_ = 0.;
0254     float sigmaRRTot_ = 0.;
0255     float varRR_ = 0.;
0256     float varZZ_ = 0.;
0257     float varEtaEta_ = 0.;
0258     float varPhiPhi_ = 0.;
0259     float sigmaRRMean_ = 0.;
0260     float sigmaZZ_ = 0.;
0261     float zBarycenter_ = 0.;
0262     float layer10percent_ = 0.;
0263     float layer50percent_ = 0.;
0264     float layer90percent_ = 0.;
0265     float triggerCells67percent_ = 0.;
0266     float triggerCells90percent_ = 0.;
0267     float first1layers_ = 0.;
0268     float first3layers_ = 0.;
0269     float first5layers_ = 0.;
0270     float firstHcal1layers_ = 0.;
0271     float firstHcal3layers_ = 0.;
0272     float firstHcal5layers_ = 0.;
0273     float last1layers_ = 0.;
0274     float last3layers_ = 0.;
0275     float last5layers_ = 0.;
0276     float emax1layers_ = 0.;
0277     float emax3layers_ = 0.;
0278     float emax5layers_ = 0.;
0279     float eot_ = 0.;
0280     int ebm0_ = 0;
0281     int ebm1_ = 0;
0282     int hbm_ = 0;
0283
0284     void updateP4AndPosition(const edm::Ptr<C>& c, bool updateCentre = true, float fraction = 1.) {
0285       double cMipt = c->mipPt() * fraction;
0286       double cPt = c->pt() * fraction;
0287       /* update cluster positions (IF requested) */
0288       if (updateCentre) {
0289         Basic3DVector<float> constituentCentre(c->position());
0290         Basic3DVector<float> clusterCentre(centre_);
0291
0292         clusterCentre = clusterCentre * mipPt_ + constituentCentre * cMipt;
0293         if ((mipPt_ + cMipt) > 0) {
0294           clusterCentre /= (mipPt_ + cMipt);
0295         }
0296         centre_ = GlobalPoint(clusterCentre);
0297
0298         if (clusterCentre.z() != 0) {
0299           centreProj_ = GlobalPoint(clusterCentre / std::abs(clusterCentre.z()));
0300         }
0301       }
0302
0303       /* update cluster energies */
0304       mipPt_ += cMipt;
0305       sumPt_ += cPt;
0306       int updatedPt = hwPt() + (int)(c->hwPt() * fraction);
0307       setHwPt(updatedPt);
0308
0309       math::PtEtaPhiMLorentzVector updatedP4(p4());
0310       updatedP4 += (c->p4() * fraction);
0311       setP4(updatedP4);
0312     }
0313   };
0314
0315 }  // namespace l1t
0316
0317 #endif