CaloTower

Macros

Line Code
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262
#ifndef DATAFORMATS_CALOTOWERS_CALOTOWER_H
#define DATAFORMATS_CALOTOWERS_CALOTOWER_H 1

#include "DataFormats/Candidate/interface/LeafCandidate.h"
#include "DataFormats/DetId/interface/DetId.h"
#include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
#include "DataFormats/Math/interface/Vector3D.h"
#include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
#include "Rtypes.h"
#include <vector>
#include <cmath>

#include "DataFormats/GeometryVector/interface/GlobalPoint.h"

/** \class CaloTower
    
\author J. Mans - Minnesota
*/

// Original author: J. Mans - Minnesota
//
// Modified: Anton Anastassov (Northwestern)
//    Make CaloTower inherit from LeafCandidate,
//    add new members and accessors.

class CaloTower : public reco::LeafCandidate {
public:
  typedef CaloTowerDetId key_type;  // for SortedCollection

  // Default constructor
  CaloTower();

  // Constructors from values

  CaloTower(const CaloTowerDetId& id,
            double emE,
            double hadE,
            double outerE,
            int ecal_tp,
            int hcal_tp,
            const PolarLorentzVector& p4,
            const GlobalPoint& emPosition,
            const GlobalPoint& hadPosition);

  CaloTower(const CaloTowerDetId& id,
            double emE,
            double hadE,
            double outerE,
            int ecal_tp,
            int hcal_tp,
            const LorentzVector& p4,
            const GlobalPoint& emPosition,
            const GlobalPoint& hadPosition);

  CaloTower(CaloTowerDetId id,
            float emE,
            float hadE,
            float outerE,
            int ecal_tp,
            int hcal_tp,
            GlobalVector p3,
            float iEnergy,
            bool massless,
            GlobalPoint emPosition,
            GlobalPoint hadPosition);

  CaloTower(CaloTowerDetId id,
            float emE,
            float hadE,
            float outerE,
            int ecal_tp,
            int hcal_tp,
            GlobalVector p3,
            float iEnergy,
            float imass,
            GlobalPoint emPosition,
            GlobalPoint hadPosition);

  // setters
  void addConstituent(DetId id) { constituents_.push_back(id); }
  void addConstituents(const std::vector<DetId>& ids);
  void setConstituents(std::vector<DetId>&& ids) { constituents_ = std::move(ids); }
  void setEcalTime(int t) { ecalTime_ = t; };
  void setHcalTime(int t) { hcalTime_ = t; };
  void setHcalSubdet(int lastHB, int lastHE, int lastHF, int lastHO) {
    int ct_ieta = ietaAbs();
    if (ct_ieta <= lastHB)
      subdet_ = HcalBarrel;
    else if (ct_ieta <= lastHE)
      subdet_ = HcalEndcap;
    else if (ct_ieta <= lastHF)
      subdet_ = HcalForward;

    //account for HO separately
    if (ct_ieta <= lastHO)
      inHO_ = true;
    else
      inHO_ = false;

    //account for gap/crossover tower separately
    if (ct_ieta == lastHB)
      inHBHEgap_ = true;
    else
      inHBHEgap_ = false;
  }

  // set CaloTower status based on the number of
  // bad/recovered/problematic cells in ECAL and HCAL

  void setCaloTowerStatus(unsigned int numBadHcalChan,
                          unsigned int numBadEcalChan,
                          unsigned int numRecHcalChan,
                          unsigned int numRecEcalChan,
                          unsigned int numProbHcalChan,
                          unsigned int numProbEcalChan);

  void setCaloTowerStatus(uint32_t s) { twrStatusWord_ = s; }

  // set the hottest cell energy in the tower
  void setHottestCellE(double e) { hottestCellE_ = e; }

  // getters
  CaloTowerDetId id() const { return id_; }
  const std::vector<DetId>& constituents() const { return constituents_; }
  size_t constituentsSize() const { return constituents_.size(); }
  DetId constituent(size_t i) const { return constituents_[i]; }

  // energy contributions from different detectors
  // energy in HO ("outerEnergy")is not included in "hadEnergy"
  double emEnergy() const { return emE_; }
  double hadEnergy() const { return hadE_; }
  double outerEnergy() const { return (inHO_) ? outerE_ : 0.0; }

  // transverse energies wrt to vtx (0,0,0)
  double emEt() const { return emE_ * sin(theta()); }
  double hadEt() const { return hadE_ * sin(theta()); }
  double outerEt() const { return (inHO_) ? outerE_ * sin(theta()) : 0.0; }

  // preserve the inherited default accessors where applicable
  // (user gets default p4 wrt to vtx (0,0,0) using p4(), etc.

  using LeafCandidate::et;
  using LeafCandidate::p;
  using LeafCandidate::p4;

  // recalculated wrt user provided vertex Z position;

  math::PtEtaPhiMLorentzVector p4(double vtxZ) const;
  double p(double vtxZ) const { return p4(vtxZ).P(); }
  double et(double vtxZ) const { return p4(vtxZ).Et(); }

  double emEt(double vtxZ) const { return emE_ * sin(p4(vtxZ).theta()); }
  double hadEt(double vtxZ) const { return hadE_ * sin(p4(vtxZ).theta()); }
  double outerEt(double vtxZ) const { return (inHO_) ? outerE_ * sin(p4(vtxZ).theta()) : 0.0; }

  // recalculated wrt vertex provided as 3D point

  math::PtEtaPhiMLorentzVector p4(const Point& v) const;
  double p(const Point& v) const { return p4(v).P(); }
  double et(const Point& v) const { return p4(v).Et(); }

  double emEt(const Point& v) const { return emE_ * sin(p4(v).theta()); }
  double hadEt(const Point& v) const { return hadE_ * sin(p4(v).theta()); }
  double outerEt(const Point& v) const { return (inHO_) ? outerE_ * sin(p4(v).theta()) : 0.0; }

  double hottestCellE() const { return hottestCellE_; }

  // Access to p4 comming from HO alone: requested by JetMET to add/subtract HO contributions
  // to the tower for cases when the tower collection was created without/with HO

  math::PtEtaPhiMLorentzVector p4_HO() const;
  math::PtEtaPhiMLorentzVector p4_HO(double vtxZ) const;
  math::PtEtaPhiMLorentzVector p4_HO(const Point& v) const;

  // the reference poins in ECAL and HCAL for direction determination
  // algorithm and parameters for selecting these points are set in the CaloTowersCreator
  const GlobalPoint& emPosition() const { return emPosition_; }
  const GlobalPoint& hadPosition() const { return hadPosition_; }

  int emLvl1() const { return emLvl1_; }
  int hadLv11() const { return hadLvl1_; }

  // energy contained in depths>1 in the HE for 18<|iEta|<29
  double hadEnergyHeOuterLayer() const { return (subdet_ == HcalEndcap) ? outerE_ : 0; }
  double hadEnergyHeInnerLayer() const { return (subdet_ == HcalEndcap) ? hadE_ - outerE_ : 0; }

  // energy in the tower by HCAL subdetector
  // This is trivial except for tower 16
  // needed by JetMET cleanup in AOD.
  double energyInHB() const;  // { return (id_.ietaAbs()<16)? outerE_ : 0.0; }
  double energyInHE() const;
  double energyInHF() const;
  double energyInHO() const;

  // time (ns) in ECAL/HCAL components of the tower based on weigted sum of the times in the contributing RecHits
  float ecalTime() const { return float(ecalTime_) * 0.01; }
  float hcalTime() const { return float(hcalTime_) * 0.01; }

  // position information on the tower
  int ieta() const { return id_.ieta(); }
  int ietaAbs() const { return id_.ietaAbs(); }
  int iphi() const { return id_.iphi(); }
  int zside() const { return id_.zside(); }

  int numCrystals() const;

  // methods to retrieve status information from the CaloTower:
  // number of bad/recovered/problematic cells in the tower
  // separately for ECAL and HCAL

  unsigned int numBadEcalCells() const { return (twrStatusWord_ & 0x1F); }
  unsigned int numRecoveredEcalCells() const { return ((twrStatusWord_ >> 5) & 0x1F); }
  unsigned int numProblematicEcalCells() const { return ((twrStatusWord_ >> 10) & 0x1F); }

  unsigned int numBadHcalCells() const { return ((twrStatusWord_ >> 15) & 0x7); }
  unsigned int numRecoveredHcalCells() const { return ((twrStatusWord_ >> 18) & 0x7); }
  unsigned int numProblematicHcalCells() const { return ((twrStatusWord_ >> 21) & 0x7); }

  // the status word itself
  uint32_t towerStatusWord() const { return twrStatusWord_; }

private:
  CaloTowerDetId id_;

  uint32_t twrStatusWord_;

  // positions of assumed EM and HAD shower positions
  GlobalPoint emPosition_;
  GlobalPoint hadPosition_;

  //hcal subdetector info
  HcalSubdetector subdet_{HcalEmpty};
  bool inHO_{false}, inHBHEgap_{false};

  // time
  int ecalTime_;
  int hcalTime_;

  float emE_, hadE_, outerE_;
  // for Jet ID use: hottest cell (ECAl or HCAL)
  float hottestCellE_;

  int emLvl1_, hadLvl1_;
  std::vector<DetId> constituents_;

  // vertex correction of EM and HAD momentum components:
  // internally used in the transformation of the CaloTower p4

  // for 3D vertex
  math::PtEtaPhiMLorentzVector hadP4(const Point& v) const;
  math::PtEtaPhiMLorentzVector emP4(const Point& v) const;

  // taking only z-component
  math::PtEtaPhiMLorentzVector hadP4(double vtxZ) const;
  math::PtEtaPhiMLorentzVector emP4(double vtxZ) const;
};

std::ostream& operator<<(std::ostream& s, const CaloTower& ct);

inline bool operator==(const CaloTower& t1, const CaloTower& t2) { return t1.id() == t2.id(); }

#endif