CaloTowers/src/CaloTower.cc

0001 #include "DataFormats/CaloTowers/interface/CaloTower.h"
0002
0003 CaloTower::CaloTower() {
0004   emE_ = 0;
0005   hadE_ = 0;
0006   outerE_ = 0;
0007   emLvl1_ = 0;
0008   hadLvl1_ = 0;
0009 }
0010
0011 CaloTower::CaloTower(const CaloTowerDetId& id,
0012                      double emE,
0013                      double hadE,
0014                      double outerE,
0015                      int ecal_tp,
0016                      int hcal_tp,
0017                      const PolarLorentzVector& p4,
0018                      const GlobalPoint& emPos,
0019                      const GlobalPoint& hadPos)
0020     : LeafCandidate(0, p4, Point(0, 0, 0)),
0021       id_(id),
0022       emPosition_(emPos),
0023       hadPosition_(hadPos),
0024       emE_(emE),
0025       hadE_(hadE),
0026       outerE_(outerE),
0027       emLvl1_(ecal_tp),
0028       hadLvl1_(hcal_tp) {}
0029
0030 CaloTower::CaloTower(const CaloTowerDetId& id,
0031                      double emE,
0032                      double hadE,
0033                      double outerE,
0034                      int ecal_tp,
0035                      int hcal_tp,
0036                      const LorentzVector& p4,
0037                      const GlobalPoint& emPos,
0038                      const GlobalPoint& hadPos)
0039     : LeafCandidate(0, p4, Point(0, 0, 0)),
0040       id_(id),
0041       emPosition_(emPos),
0042       hadPosition_(hadPos),
0043       emE_(emE),
0044       hadE_(hadE),
0045       outerE_(outerE),
0046       emLvl1_(ecal_tp),
0047       hadLvl1_(hcal_tp) {}
0048
0049 CaloTower::CaloTower(CaloTowerDetId id,
0050                      float emE,
0051                      float hadE,
0052                      float outerE,
0053                      int ecal_tp,
0054                      int hcal_tp,
0055                      GlobalVector p3,
0056                      float iEnergy,
0057                      bool massless,
0058                      GlobalPoint emPos,
0059                      GlobalPoint hadPos)
0060     : LeafCandidate(0, p3, iEnergy, massless, Point(0, 0, 0)),
0061       id_(id),
0062       emPosition_(emPos),
0063       hadPosition_(hadPos),
0064       emE_(emE),
0065       hadE_(hadE),
0066       outerE_(outerE),
0067       emLvl1_(ecal_tp),
0068       hadLvl1_(hcal_tp) {}
0069
0070 CaloTower::CaloTower(CaloTowerDetId id,
0071                      float emE,
0072                      float hadE,
0073                      float outerE,
0074                      int ecal_tp,
0075                      int hcal_tp,
0076                      GlobalVector p3,
0077                      float iEnergy,
0078                      float imass,
0079                      GlobalPoint emPos,
0080                      GlobalPoint hadPos)
0081     : LeafCandidate(0, p3, iEnergy, imass, Point(0, 0, 0)),
0082       id_(id),
0083       emPosition_(emPos),
0084       hadPosition_(hadPos),
0085       emE_(emE),
0086       hadE_(hadE),
0087       outerE_(outerE),
0088       emLvl1_(ecal_tp),
0089       hadLvl1_(hcal_tp) {}
0090
0091 // recalculated momentum-related quantities wrt user provided vertex Z position
0092
0093 math::PtEtaPhiMLorentzVector CaloTower::hadP4(double vtxZ) const {
0094   // note: for now we use the same position for HO as for the other detectors
0095
0096   double hcalTot;
0097   if (inHO_)
0098     hcalTot = (energy() - emE_);
0099   else
0100     hcalTot = hadE_;
0101
0102   if (hcalTot > 0) {
0103     double ctgTheta = (hadPosition_.z() - vtxZ) / hadPosition_.perp();
0104     double newEta = asinh(ctgTheta);
0105     double pf = 1.0 / cosh(newEta);
0106
0107     return PolarLorentzVector(hcalTot * pf, newEta, hadPosition_.phi(), 0.0);
0108   }
0109
0110   return math::PtEtaPhiMLorentzVector(0, 0, 0, 0);
0111 }
0112
0113 math::PtEtaPhiMLorentzVector CaloTower::emP4(double vtxZ) const {
0114   if (emE_ > 0) {
0115     double ctgTheta = (emPosition_.z() - vtxZ) / emPosition_.perp();
0116     double newEta = asinh(ctgTheta);
0117     double pf = 1.0 / cosh(newEta);
0118
0119     return math::PtEtaPhiMLorentzVector(emE_ * pf, newEta, emPosition_.phi(), 0.0);
0120   }
0121
0122   return math::PtEtaPhiMLorentzVector(0, 0, 0, 0);
0123 }
0124
0125 // recalculated momentum-related quantities wrt user provided 3D vertex
0126
0127 math::PtEtaPhiMLorentzVector CaloTower::hadP4(const Point& v) const {
0128   // note: for now we use the same position for HO as for the other detectors
0129
0130   double hcalTot;
0131   if (inHO_)
0132     hcalTot = (energy() - emE_);
0133   else
0134     hcalTot = hadE_;
0135
0136   if (hcalTot > 0) {
0137     GlobalPoint p(v.x(), v.y(), v.z());
0138     math::XYZVector dir = math::XYZVector(hadPosition_ - p);
0139     return math::PtEtaPhiMLorentzVector(hcalTot * sin(dir.theta()), dir.eta(), dir.phi(), 0.0);
0140   }
0141
0142   return math::PtEtaPhiMLorentzVector(0, 0, 0, 0);
0143 }
0144
0145 math::PtEtaPhiMLorentzVector CaloTower::emP4(const Point& v) const {
0146   if (emE_ > 0) {
0147     GlobalPoint p(v.x(), v.y(), v.z());
0148     math::XYZVector dir = math::XYZVector(emPosition_ - p);
0149     return math::PtEtaPhiMLorentzVector(emE_ * sin(dir.theta()), dir.eta(), dir.phi(), 0.0);
0150   }
0151
0152   return math::PtEtaPhiMLorentzVector(0, 0, 0, 0);
0153 }
0154
0155 math::PtEtaPhiMLorentzVector CaloTower::p4(double vtxZ) const {
0156   if (subdet_ == HcalBarrel || subdet_ == HcalEndcap) {
0157     return (emP4(vtxZ) + hadP4(vtxZ));
0158   }
0159   // em and had energy in HF are defined in a special way
0160   double ctgTheta =
0161       (emPosition_.z() - vtxZ) / emPosition_.perp();  // em and had positions in HF are forced to be the same
0162   double newEta = asinh(ctgTheta);
0163   double pf = 1.0 / cosh(newEta);
0164   return math::PtEtaPhiMLorentzVector(p4().energy() * pf, newEta, emPosition_.phi(), 0.0);
0165 }
0166
0167 math::PtEtaPhiMLorentzVector CaloTower::p4(const Point& v) const {
0168   if (subdet_ == HcalBarrel || subdet_ == HcalEndcap) {
0169     return emP4(v) + hadP4(v);
0170   }
0171   // em and had energy in HF are defined in a special way
0172   GlobalPoint p(v.x(), v.y(), v.z());
0173   math::XYZVector dir = math::XYZVector(emPosition_ - p);  // em and had positions in HF are forced to be the same
0174   return math::PtEtaPhiMLorentzVector(p4().energy() * sin(dir.theta()), dir.eta(), dir.phi(), 0.0);
0175 }
0176
0177 // p4 contribution from HO alone (note: direction is always taken to be the same as used for HB.)
0178
0179 math::PtEtaPhiMLorentzVector CaloTower::p4_HO(const Point& v) const {
0180   if (!inHO_ || outerE_ < 0)
0181     return math::PtEtaPhiMLorentzVector(0, 0, 0, 0);
0182
0183   GlobalPoint p(v.x(), v.y(), v.z());
0184   math::XYZVector dir = math::XYZVector(hadPosition_ - p);
0185   return math::PtEtaPhiMLorentzVector(outerE_ * sin(dir.theta()), dir.eta(), dir.phi(), 0.0);
0186 }
0187
0188 math::PtEtaPhiMLorentzVector CaloTower::p4_HO(double vtxZ) const {
0189   Point p(0, 0, vtxZ);
0190   return p4_HO(p);
0191 }
0192
0193 math::PtEtaPhiMLorentzVector CaloTower::p4_HO() const {
0194   if (!inHO_ || outerE_ < 0)
0195     return math::PtEtaPhiMLorentzVector(0.0, 0.0, 0.0, 0.0);
0196   return math::PtEtaPhiMLorentzVector(outerE_ * sin(hadPosition_.theta()), hadPosition_.eta(), hadPosition_.phi(), 0.0);
0197 }
0198
0199 void CaloTower::addConstituents(const std::vector<DetId>& ids) {
0200   constituents_.reserve(constituents_.size() + ids.size());
0201   constituents_.insert(constituents_.end(), ids.begin(), ids.end());
0202 }
0203
0204 int CaloTower::numCrystals() const {
0205   if (subdet_ == HcalForward)
0206     return 0;
0207
0208   int nC = 0;
0209   std::vector<DetId>::const_iterator it = constituents_.begin();
0210   for (; it != constituents_.end(); ++it) {
0211     if (it->det() == DetId::Ecal)
0212       ++nC;
0213   }
0214
0215   return nC;
0216 }
0217
0218 // Set the CaloTower status word from the number of bad/recovered/problematic
0219 // cells in HCAL and ECAL.
0220
0221 void CaloTower::setCaloTowerStatus(unsigned int numBadHcalChan,
0222                                    unsigned int numBadEcalChan,
0223                                    unsigned int numRecHcalChan,
0224                                    unsigned int numRecEcalChan,
0225                                    unsigned int numProbHcalChan,
0226                                    unsigned int numProbEcalChan) {
0227   twrStatusWord_ = 0x0;
0228
0229   twrStatusWord_ |= (numBadEcalChan & 0x1F);
0230   twrStatusWord_ |= ((numRecEcalChan & 0x1F) << 5);
0231   twrStatusWord_ |= ((numProbEcalChan & 0x1F) << 10);
0232   twrStatusWord_ |= ((numBadHcalChan & 0x7) << 15);
0233   twrStatusWord_ |= ((numRecHcalChan & 0x7) << 18);
0234   twrStatusWord_ |= ((numProbHcalChan & 0x7) << 21);
0235
0236   return;
0237 }
0238
0239 // energy in the tower by HCAL subdetector
0240 // This is trivia except for tower 16
0241 // needed by JetMET cleanup in AOD.
0242
0243 double CaloTower::energyInHB() const {
0244   if (inHBHEgap_)
0245     return hadE_ - outerE_;
0246   else if (subdet_ == HcalBarrel)
0247     return hadE_;
0248   else
0249     return 0.0;
0250 }
0251
0252 double CaloTower::energyInHE() const {
0253   if (inHBHEgap_)
0254     return outerE_;
0255   else if (subdet_ == HcalEndcap)
0256     return hadE_;
0257   else
0258     return 0.0;
0259 }
0260
0261 double CaloTower::energyInHF() const {
0262   if (subdet_ == HcalForward)
0263     return energy();
0264   else
0265     return 0.0;
0266 }
0267
0268 // this is actual energy contributed to the tower
0269 // (outerEnergy() returns HO energy regardless if it is used or not)
0270 // Note: rounding error may lead to values not identically equal to zero
0271 // when HO was not used
0272
0273 double CaloTower::energyInHO() const {
0274   if (inHO_)
0275     return (energy() - hadE_ - emE_);
0276   else
0277     return 0.0;
0278 }
0279
0280 std::ostream& operator<<(std::ostream& s, const CaloTower& ct) {
0281   return s << ct.id() << ":" << ct.et() << " GeV ET (EM=" << ct.emEt() << " HAD=" << ct.hadEt()
0282            << " OUTER=" << ct.outerEt() << ") (" << ct.eta() << "," << ct.phi() << ")";
0283 }