EcalRecHit/interface/EcalRecHit.h

0001 #ifndef DATAFORMATS_ECALRECHIT_H
0002 #define DATAFORMATS_ECALRECHIT_H 1
0003
0004 #include "DataFormats/DetId/interface/DetId.h"
0005 #include "DataFormats/CaloRecHit/interface/CaloRecHit.h"
0006
0007 #include <vector>
0008 #include <cmath>
0009
0010 /** \class EcalRecHit
0011  *
0012  * \author P. Meridiani INFN Roma1
0013  */
0014
0015 class EcalRecHit {
0016 public:
0017   typedef DetId key_type;
0018
0019   // recHit flags
0020   enum Flags {
0021     kGood = 0,   // channel ok, the energy and time measurement are reliable
0022     kPoorReco,   // the energy is available from the UncalibRecHit, but approximate (bad shape, large chi2)
0023     kOutOfTime,  // the energy is available from the UncalibRecHit (sync reco), but the event is out of time
0024     kFaultyHardware,  // The energy is available from the UncalibRecHit, channel is faulty at some hardware level (e.g. noisy)
0025     kNoisy,           // the channel is very noisy
0026     kPoorCalib,  // the energy is available from the UncalibRecHit, but the calibration of the channel is poor
0027     kSaturated,  // saturated channel (recovery not tried)
0028     kLeadingEdgeRecovered,  // saturated channel: energy estimated from the leading edge before saturation
0029     kNeighboursRecovered,   // saturated/isolated dead: energy estimated from neighbours
0030     kTowerRecovered,        // channel in TT with no data link, info retrieved from Trigger Primitive
0031     kDead,                  // channel is dead and any recovery fails
0032     kKilled,                // MC only flag: the channel is killed in the real detector
0033     kTPSaturated,           // the channel is in a region with saturated TP
0034     kL1SpikeFlag,           // the channel is in a region with TP with sFGVB = 0
0035     kWeird,                 // the signal is believed to originate from an anomalous deposit (spike)
0036     kDiWeird,               // the signal is anomalous, and neighbors another anomalous signal
0037     kHasSwitchToGain6,      // at least one data frame is in G6
0038     kHasSwitchToGain1,      // at least one data frame is in G1
0039                             //
0040     kUnknown                // to ease the interface with functions returning flags.
0041   };
0042
0043   // ES recHit flags
0044   enum ESFlags {
0045     kESGood,
0046     kESDead,
0047     kESHot,
0048     kESPassBX,
0049     kESTwoGoodRatios,
0050     kESBadRatioFor12,
0051     kESBadRatioFor23Upper,
0052     kESBadRatioFor23Lower,
0053     kESTS1Largest,
0054     kESTS3Largest,
0055     kESTS3Negative,
0056     kESSaturated,
0057     kESTS2Saturated,
0058     kESTS3Saturated,
0059     kESTS13Sigmas,
0060     kESTS15Sigmas
0061   };
0062
0063   EcalRecHit() : energy_(0), time_(0), flagBits_(0) {}
0064   // by default a recHit is greated with no flag
0065   explicit EcalRecHit(const DetId& id, float energy, float time, uint32_t extra = 0, uint32_t flagBits = 0)
0066       : id_(id), energy_(energy), time_(time), flagBits_(flagBits), extra_(extra) {}
0067
0068   float energy() const { return energy_; }
0069   void setEnergy(float energy) { energy_ = energy; }
0070   float time() const { return time_; }
0071   void setTime(float time) { time_ = time; }
0072   const DetId& detid() const { return id_; }
0073
0074   /// get the id
0075   // For the moment not returning a specific id for subdetector
0076   // @TODO why this method?! should use detid()
0077   DetId id() const { return DetId(detid()); }
0078
0079   bool isRecovered() const {
0080     return checkFlag(kLeadingEdgeRecovered) || checkFlag(kNeighboursRecovered) || checkFlag(kTowerRecovered);
0081   }
0082
0083   bool isTimeValid() const { return (this->timeError() > 0); }
0084
0085   bool isTimeErrorValid() const {
0086     if (!isTimeValid())
0087       return false;
0088
0089     if (timeError() >= 10000)
0090       return false;
0091
0092     return true;
0093   }
0094
0095   static inline uint32_t getMasked(uint32_t value, uint32_t offset, uint32_t width) {
0096     return (value >> offset) & ((1 << width) - 1);
0097   }
0098
0099   static inline uint32_t setMasked(uint32_t value, uint32_t x, uint32_t offset, uint32_t width) {
0100     const uint32_t mask = ((1 << width) - 1) << offset;
0101     value &= ~mask;
0102     value |= (x & ((1U << width) - 1)) << offset;
0103     return value;
0104   }
0105
0106   static inline int getPower10(float e) {
0107     static constexpr float p10[] = {1.e-2f, 1.e-1f, 1.f, 1.e1f, 1.e2f, 1.e3f, 1.e4f, 1.e5f, 1.e6f};
0108     int b = e < p10[4] ? 0 : 5;
0109     for (; b < 9; ++b)
0110       if (e < p10[b])
0111         break;
0112     return b;
0113   }
0114
0115   /* the new bit structure
0116    * 0..6   - chi2 in time events (chi2()), offset=0, width=7
0117    * 8..20  - energy uncertainty, offset=8, width=13
0118    * 24..31 - timeError(), offset=24, width=8
0119    */
0120   float chi2() const {
0121     uint32_t rawChi2 = getMasked(extra_, 0, 7);
0122     return (float)rawChi2 / (float)((1 << 7) - 1) * 64.;
0123   }
0124
0125   void setChi2(float chi2) {
0126     // bound the max value of the chi2
0127     if (chi2 > 64)
0128       chi2 = 64;
0129
0130     // use 7 bits
0131     uint32_t rawChi2 = lround(chi2 / 64. * ((1 << 7) - 1));
0132     extra_ = setMasked(extra_, rawChi2, 0, 7);
0133   }
0134
0135   float energyError() const {
0136     uint32_t rawEnergy = getMasked(extra_, 8, 13);
0137     uint16_t exponent = rawEnergy >> 10;
0138     uint16_t significand = ~(0xE << 9) & rawEnergy;
0139     return (float)significand * pow(10, exponent - 5);
0140   }
0141
0142   // set the energy uncertainty
0143   // (only energy >= 0 will be stored)
0144   void setEnergyError(float energy) {
0145     uint32_t rawEnergy = 0;
0146     if (energy > 0.001) {
0147       uint16_t exponent = getPower10(energy);
0148       static constexpr float ip10[] = {1.e5f, 1.e4f, 1.e3f, 1.e2f, 1.e1f, 1.e0f, 1.e-1f, 1.e-2f, 1.e-3f, 1.e-4};
0149       uint16_t significand = lround(energy * ip10[exponent]);
0150       // use 13 bits (3 exponent, 10 significand)
0151       rawEnergy = exponent << 10 | significand;
0152       /* here for doc and regression test
0153       uint16_t exponent_old = lround(floor(log10(energy))) + 3;
0154       uint16_t significand_old = lround(energy/pow(10, exponent - 5));
0155       std::cout << energy << ' ' << exponent << ' ' << significand
0156                           << ' ' << exponent_old << ' ' << significand_old << std::endl;
0157       assert(exponent==exponent_old);
0158       assert(significand==significand_old);
0159       */
0160     }
0161
0162     extra_ = setMasked(extra_, rawEnergy, 8, 13);
0163   }
0164
0165   float timeError() const {
0166     uint32_t timeErrorBits = getMasked(extra_, 24, 8);
0167     // all bits off --> time reco bailed out (return negative value)
0168     if ((0xFF & timeErrorBits) == 0x00)
0169       return -1;
0170     // all bits on --> time error over 5 ns (return large value)
0171     if ((0xFF & timeErrorBits) == 0xFF)
0172       return 10000;
0173
0174     float LSB = 1.26008;
0175     uint8_t exponent = timeErrorBits >> 5;
0176     uint8_t significand = timeErrorBits & ~(0x7 << 5);
0177     return pow(2., exponent) * significand * LSB / 1000.;
0178   }
0179
0180   void setTimeError(uint8_t timeErrBits) { extra_ = setMasked(extra_, timeErrBits & 0xFF, 24, 8); }
0181
0182   /// set the flags (from Flags or ESFlags)
0183   void setFlag(int flag) { flagBits_ |= (0x1 << flag); }
0184   void unsetFlag(int flag) { flagBits_ &= ~(0x1 << flag); }
0185
0186   /// check if the flag is true
0187   bool checkFlag(int flag) const { return flagBits_ & (0x1 << flag); }
0188
0189   /// check if one of the flags in a set is true
0190   bool checkFlags(const std::vector<int>& flagsvec) const {
0191     for (std::vector<int>::const_iterator flagPtr = flagsvec.begin(); flagPtr != flagsvec.end();
0192          ++flagPtr) {  // check if one of the flags is up
0193
0194       if (checkFlag(*flagPtr))
0195         return true;
0196     }
0197
0198     return false;
0199   }
0200
0201   uint32_t flagsBits() const { return flagBits_; }
0202
0203   /// apply a bitmask to our flags. Experts only
0204   bool checkFlagMask(uint32_t mask) const { return flagBits_ & mask; }
0205
0206   /// DEPRECATED provided for temporary backward compatibility
0207   Flags recoFlag() const {
0208     for (int i = kUnknown;; --i) {
0209       if (checkFlag(i))
0210         return Flags(i);
0211       if (i == 0)
0212         break;
0213     }
0214
0215     // no flag assigned, assume good
0216     return kGood;
0217   }
0218
0219 private:
0220   // from calorechit
0221   DetId id_;
0222   float energy_;
0223   float time_;
0224
0225   /// store rechit condition (see Flags enum) in a bit-wise way
0226   uint32_t flagBits_;
0227
0228   // packed uint32_t for timeError, chi2, energyError
0229   uint32_t extra_;
0230 };
0231
0232 std::ostream& operator<<(std::ostream& s, const EcalRecHit& hit);
0233
0234 #endif