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