HGCRecHit/src/HGCRecHit.cc

0001 #include "DataFormats/HGCRecHit/interface/HGCRecHit.h"
0002 #include "DataFormats/ForwardDetId/interface/HFNoseDetId.h"
0003 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0004 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0005 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0006 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0007 #include <cassert>
0008 #include <cmath>
0009
0010 HGCRecHit::HGCRecHit() : CaloRecHit(), flagBits_(0) {}
0011
0012 HGCRecHit::HGCRecHit(
0013     const DetId& id, float energy, float time, uint32_t flags, uint32_t flagBits, uint8_t son, float timeError)
0014     : CaloRecHit(id, energy, time, flags), flagBits_(flagBits), signalOverSigmaNoise_(son), timeError_(timeError) {}
0015
0016 float HGCRecHit::chi2() const {
0017   uint32_t rawChi2 = 0x7F & (flags() >> 4);
0018   return (float)rawChi2 / (float)((1 << 7) - 1) * 64.f;
0019 }
0020
0021 float HGCRecHit::outOfTimeChi2() const {
0022   uint32_t rawChi2Prob = 0x7F & (flags() >> 24);
0023   return (float)rawChi2Prob / (float)((1 << 7) - 1) * 64.f;
0024 }
0025
0026 float HGCRecHit::outOfTimeEnergy() const {
0027   uint32_t rawEnergy = (0x1FFF & flags() >> 11);
0028   uint16_t exponent = rawEnergy >> 10;
0029   uint16_t significand = ~(0xE << 9) & rawEnergy;
0030   return (float)significand * pow(10, exponent - 5);
0031 }
0032
0033 void HGCRecHit::setChi2(float chi2) {
0034   // bound the max value of the chi2
0035   if (chi2 > 64)
0036     chi2 = 64;
0037   // use 7 bits
0038   uint32_t rawChi2 = lround(chi2 / 64.f * ((1 << 7) - 1));
0039   // shift by 4 bits (recoFlag)
0040   setFlags((~(0x7F << 4) & flags()) | ((rawChi2 & 0x7F) << 4));
0041 }
0042
0043 void HGCRecHit::setOutOfTimeEnergy(float energy) {
0044   if (energy > 0.001f) {
0045     uint16_t exponent = lround(floor(log10(energy))) + 3;
0046     uint16_t significand = lround(energy / pow(10, exponent - 5));
0047     // use 13 bits (3 exponent, 10 significand)
0048     uint32_t rawEnergy = exponent << 10 | significand;
0049     // shift by 11 bits (recoFlag + chi2)
0050     setFlags((~(0x1FFF << 11) & flags()) | ((rawEnergy & 0x1FFF) << 11));
0051   }
0052 }
0053
0054 void HGCRecHit::setOutOfTimeChi2(float chi2) {
0055   // bound the max value of chi2
0056   if (chi2 > 64)
0057     chi2 = 64;
0058   // use 7 bits
0059   uint32_t rawChi2 = lround(chi2 / 64.f * ((1 << 7) - 1));
0060   // shift by 24 bits (recoFlag + chi2 + outOfTimeEnergy)
0061   setFlags((~(0x7F << 24) & flags()) | ((rawChi2 & 0x7F) << 24));
0062 }
0063
0064 void HGCRecHit::setSignalOverSigmaNoise(float sOverNoise) {
0065   // bound the max value of sOverNoise
0066   if (sOverNoise > 32.f)
0067     sOverNoise = 32.f;
0068   //use 8 bits
0069   signalOverSigmaNoise_ = lround(sOverNoise / 32.f * ((1 << 8) - 1));
0070 }
0071
0072 float HGCRecHit::signalOverSigmaNoise() const { return (float)signalOverSigmaNoise_ * 0.125f; }
0073
0074 void HGCRecHit::setTimeError(float timeErr) {
0075   //expected resolution on single cell for that given S/N
0076   timeError_ = timeErr;
0077 }
0078
0079 float HGCRecHit::timeError() const { return timeError_; }
0080
0081 bool HGCRecHit::isTimeValid() const {
0082   if (timeError() <= 0)
0083     return false;
0084   else
0085     return true;
0086 }
0087
0088 bool HGCRecHit::isTimeErrorValid() const {
0089   if (!isTimeValid())
0090     return false;
0091   if (timeError() >= 10000)
0092     return false;
0093
0094   return true;
0095 }
0096
0097 /// check if one of the flags in a set is true
0098 bool HGCRecHit::checkFlags(const std::vector<int>& flagsvec) const {
0099   for (std::vector<int>::const_iterator flagPtr = flagsvec.begin(); flagPtr != flagsvec.end();
0100        ++flagPtr) {  // check if one of the flags is up
0101     if (checkFlag(*flagPtr))
0102       return true;
0103   }
0104   return false;
0105 }
0106
0107 std::ostream& operator<<(std::ostream& s, const HGCRecHit& hit) {
0108   if (hit.detid().det() == DetId::Forward && hit.detid().subdetId() == HGCEE)
0109     return s << HGCalDetId(hit.detid()) << ": " << hit.energy() << " GeV, " << hit.time() << " ns";
0110   else if (hit.detid().det() == DetId::Forward && hit.detid().subdetId() == HGCHEF)
0111     return s << HGCalDetId(hit.detid()) << ": " << hit.energy() << " GeV, " << hit.time() << " ns";
0112   else if (hit.detid().det() == DetId::Hcal && hit.detid().subdetId() == HcalEndcap)
0113     return s << HcalDetId(hit.detid()) << ": " << hit.energy() << " GeV, " << hit.time() << " ns";
0114   else if (hit.detid().det() == DetId::HGCalEE || hit.detid().det() == DetId::HGCalHSi)
0115     return s << HGCSiliconDetId(hit.detid()) << ": " << hit.energy() << " GeV, " << hit.time() << " ns";
0116   else if (hit.detid().det() == DetId::HGCalHSc)
0117     return s << HGCScintillatorDetId(hit.detid()) << ": " << hit.energy() << " GeV, " << hit.time() << " ns";
0118   else if (hit.detid().det() == DetId::Forward && hit.detid().subdetId() == HFNose)
0119     return s << HFNoseDetId(hit.detid()) << ": " << hit.energy() << " GeV, " << hit.time() << " ns";
0120   else
0121     return s << "HGCRecHit undefined subdetector";
0122 }