EgammaTools/interface/EnergyScaleCorrection.h

0001 #ifndef RecoEgamma_EgammaTools_EnergyScaleCorrection_h
0002 #define RecoEgamma_EgammaTools_EnergyScaleCorrection_h
0003
0004 //author: Alan Smithee
0005 //description:
0006 //  A port of Shervin Nourbakhsh's EnergyScaleCorrection_class in EgammaAnalysis/ElectronTools
0007 //  this reads the scale & smearing corrections in from a text file for given categories
0008 //  it then allows these values to be accessed
0009
0010 #include <iostream>
0011 #include <fstream>
0012 #include <vector>
0013 #include <map>
0014 #include <cmath>
0015 #include <string>
0016 #include <bitset>
0017
0018 class EnergyScaleCorrection {
0019 public:
0020   enum FileFormat { UNKNOWN = 0, GLOBE, ECALELF_TOY, ECALELF };
0021
0022   enum ParamSmear { kNone = 0, kRho, kPhi, kNParamSmear };
0023
0024   enum ScaleNuisances {
0025     kErrStatBitNr = 0,
0026     kErrSystBitNr = 1,
0027     kErrGainBitNr = 2,
0028     kErrNrBits = 3,
0029     kErrNone = 0,
0030     kErrStat = 1,
0031     kErrSyst = 2,
0032     kErrGain = 4,
0033     kErrStatSyst = 3,
0034     kErrStatGain = 5,
0035     kErrSystGain = 6,
0036     kErrStatSystGain = 7
0037   };
0038
0039   class ScaleCorrection {
0040   public:
0041     ScaleCorrection() : scale_(1.), scaleErrStat_(0.), scaleErrSyst_(0.), scaleErrGain_(0.) {}
0042     ScaleCorrection(float iScale, float iScaleErrStat, float iScaleErrSyst, float iScaleErrGain)
0043         : scale_(iScale), scaleErrStat_(iScaleErrStat), scaleErrSyst_(iScaleErrSyst), scaleErrGain_(iScaleErrGain) {}
0044
0045     float scale() const { return scale_; }
0046     float scaleErr(const std::bitset<kErrNrBits>& uncBitMask) const;
0047     float scaleErrStat() const { return scaleErrStat_; }
0048     float scaleErrSyst() const { return scaleErrSyst_; }
0049     float scaleErrGain() const { return scaleErrGain_; }
0050
0051     friend std::ostream& operator<<(std::ostream& os, const ScaleCorrection& a) { return a.print(os); }
0052     std::ostream& print(std::ostream& os) const;
0053
0054   private:
0055     float scale_, scaleErrStat_, scaleErrSyst_, scaleErrGain_;
0056   };
0057
0058   struct SmearCorrection {
0059   public:
0060     SmearCorrection() : rho_(0.), rhoErr_(0.), phi_(0.), phiErr_(0.), eMean_(0.), eMeanErr_(0.) {}
0061     SmearCorrection(float iRho, float iRhoErr, float iPhi, float iPhiErr, float iEMean, float iEMeanErr)
0062         : rho_(iRho), rhoErr_(iRhoErr), phi_(iPhi), phiErr_(iPhiErr), eMean_(iEMean), eMeanErr_(iEMeanErr) {}
0063
0064     friend std::ostream& operator<<(std::ostream& os, const SmearCorrection& a) { return a.print(os); }
0065     std::ostream& print(std::ostream& os) const;
0066
0067     float sigma(const float et, const float nrSigmaRho = 0., const float nrSigmaPhi = 0.) const {
0068       const float rhoVal = rho_ + rhoErr_ * nrSigmaRho;
0069       const float phiVal = phi_ + phiErr_ * nrSigmaPhi;
0070       const float constTerm = rhoVal * std::sin(phiVal);
0071       const float alpha = rhoVal * eMean_ * std::cos(phiVal);
0072       return std::sqrt(constTerm * constTerm + alpha * alpha / et);
0073     }
0074
0075   private:
0076     float rho_, rhoErr_;
0077     float phi_, phiErr_;
0078     float eMean_, eMeanErr_;
0079   };
0080
0081   class CorrectionCategory {
0082   public:
0083     CorrectionCategory(const std::string& category, int runnrMin = 0, int runnrMax = 999999);
0084     CorrectionCategory(
0085         const unsigned int runnr, const float et, const float eta, const float r9, const unsigned int gainSeed)
0086         : runMin_(runnr),
0087           runMax_(runnr),
0088           etaMin_(std::abs(eta)),
0089           etaMax_(std::abs(eta)),
0090           r9Min_(r9),
0091           r9Max_(r9),
0092           etMin_(et),
0093           etMax_(et),
0094           gain_(gainSeed) {}
0095
0096     CorrectionCategory(unsigned int runMin,
0097                        unsigned int runMax,
0098                        float etaMin,
0099                        float etaMax,
0100                        float r9Min,
0101                        float r9Max,
0102                        float etMin,
0103                        float etMax,
0104                        unsigned int gainSeed);
0105
0106     bool operator<(const CorrectionCategory& b) const;
0107     bool inCategory(
0108         const unsigned int runnr, const float et, const float eta, const float r9, const unsigned int gainSeed) const;
0109
0110     friend std::ostream& operator<<(std::ostream& os, const CorrectionCategory& a) { return a.print(os); }
0111     std::ostream& print(std::ostream& os) const;
0112
0113   private:
0114     //all boundaries are inclusive (X<=Y<=Z)
0115     unsigned int runMin_;
0116     unsigned int runMax_;
0117     float etaMin_;       ///< min eta value for the bin
0118     float etaMax_;       ///< max eta value for the bin
0119     float r9Min_;        ///< min R9 vaule for the bin
0120     float r9Max_;        ///< max R9 value for the bin
0121     float etMin_;        ///< min Et value for the bin
0122     float etMax_;        ///< max Et value for the bin
0123     unsigned int gain_;  ///< 12, 6, 1, 61 (double gain switch)
0124   };
0125
0126 public:
0127   EnergyScaleCorrection(const std::string& correctionFileName, unsigned int genSeed = 0);
0128   EnergyScaleCorrection() {}
0129   ~EnergyScaleCorrection() {}
0130
0131   float scaleCorr(unsigned int runnr,
0132                   double et,
0133                   double eta,
0134                   double r9,
0135                   unsigned int gainSeed = 12,
0136                   std::bitset<kErrNrBits> uncBitMask = kErrNone) const;
0137
0138   float scaleCorrUncert(unsigned int runnr,
0139                         double et,
0140                         double eta,
0141                         double r9,
0142                         unsigned int gainSeed,
0143                         std::bitset<kErrNrBits> uncBitMask = kErrNone) const;
0144
0145   float smearingSigma(
0146       int runnr, double et, double eta, double r9, unsigned int gainSeed, ParamSmear par, float nSigma = 0.) const;
0147   float smearingSigma(
0148       int runnr, double et, double eta, double r9, unsigned int gainSeed, float nSigmaRho, float nSigmaPhi) const;
0149
0150   void setSmearingType(FileFormat value);
0151
0152   const ScaleCorrection* getScaleCorr(unsigned int runnr, double et, double eta, double r9, unsigned int gainSeed) const;
0153   const SmearCorrection* getSmearCorr(unsigned int runnr, double et, double eta, double r9, unsigned int gainSeed) const;
0154
0155 private:
0156   void addScale(const std::string& category,
0157                 int runMin,
0158                 int runMax,
0159                 double deltaP,
0160                 double errDeltaP,
0161                 double errSystDeltaP,
0162                 double errDeltaPGain);
0163
0164   void addScale(int runMin,
0165                 int runMax,
0166                 double etaMin,
0167                 double etaMax,
0168                 double r9Min,
0169                 double r9Max,
0170                 double etMin,
0171                 double etMax,
0172                 unsigned int gain,
0173                 double energyScale,
0174                 double energyScaleErrStat,
0175                 double energyScaleErrSyst,
0176                 double energyScaleErrGain);
0177
0178   void addSmearing(const std::string& category,
0179                    int runMin,
0180                    int runMax,
0181                    double rho,
0182                    double errRho,
0183                    double phi,
0184                    double errPhi,
0185                    double eMean,
0186                    double errEMean);
0187
0188   void readScalesFromFile(const std::string& filename);
0189   void readSmearingsFromFile(const std::string& filename);
0190
0191   //static data members
0192   static constexpr float kDefaultScaleVal_ = 1.0;
0193   static constexpr float kDefaultSmearVal_ = 0.0;
0194
0195   //data members
0196   FileFormat smearingType_;
0197   std::map<CorrectionCategory, ScaleCorrection> scales_;
0198   std::map<CorrectionCategory, SmearCorrection> smearings_;
0199
0200   template <typename T1, typename T2>
0201   class Sorter {
0202   public:
0203     bool operator()(const std::pair<T1, T2>& lhs, const T1& rhs) const { return lhs.first < rhs; }
0204     bool operator()(const std::pair<T1, T2>& lhs, const std::pair<T1, T2>& rhs) const { return lhs.first < rhs.first; }
0205     bool operator()(const T1& lhs, const std::pair<T1, T2>& rhs) const { return lhs < rhs.first; }
0206     bool operator()(const T1& lhs, const T1& rhs) const { return lhs < rhs; }
0207   };
0208 };
0209
0210 #endif