LightFlavorRatios/systematic_uncertainty/SystematicUncertainty.h

0001 #ifndef SYSTEMATIC_UNCERTAINTY_H
0002 #define SYSTEMATIC_UNCERTAINTY_H
0003
0004 class SystematicUncertainty
0005 {
0006   public:
0007
0008   SystematicUncertainty(TMatrixDSym& cov_in, std::vector<double>& bins_in)
0009   {
0010     covariance = cov_in;
0011     bins = bins_in;
0012     hBins = std::make_unique<TH1D>("hBins","hBins",bins.size()-1,bins.data());
0013   }
0014
0015   bool in_range(double x)
0016   {
0017     return (x>=bins[0] && x<=bins[bins.size()-1]);
0018   }
0019
0020   void check_range_warn(double x_low, double x_high)
0021   {
0022     if(!in_range(x_low) || !in_range(x_high))
0023     {
0024       std::cout << "WARNING: selection [" << x_low << ", " << x_high << "]"
0025         << " is out of range [" << get_xmin() << ", "
0026         << get_xmax() << "]" << std::endl;
0027     }
0028   }
0029
0030   double uncorrelatedUncertaintyAt(double x)
0031   {
0032     if(!in_range(x))
0033     {
0034       std::cout << "WARNING: point " << x << " is out of range [" << bins[0] << ", " << bins[bins.size()-1] << "]" << std::endl;
0035       std::cout << "returning 0" << std::endl;
0036       return 0.;
0037     }
0038
0039     int iBin = hBins->FindBin(x)-1;
0040     return sqrt(covariance(iBin,iBin));
0041   }
0042
0043   double nextBinCovarianceAt(double x)
0044   {
0045     if(!in_range(x))
0046     {
0047       std::cout << "WARNING: point " << x << " is out of range [" << bins[0] << ", " << bins[bins.size()-1] << "]" << std::endl;
0048       std::cout << "returning 0" << std::endl;
0049       return 0.;
0050     }
0051
0052     int iBin = hBins->FindBin(x)-1;
0053     if(iBin<bins.size()) return covariance(iBin,iBin+1);
0054     else return 0.;
0055   }
0056
0057   double prevBinCovarianceAt(double x)
0058   {
0059     if(!in_range(x))
0060     {
0061       std::cout << "WARNING: point " << x << " is out of range [" << bins[0] << ", " << bins[bins.size()-1] << "]" << std::endl;
0062       std::cout << "returning 0" << std::endl;
0063       return 0.;
0064     }
0065
0066     int iBin = hBins->FindBin(x)-1;
0067     if(iBin>0) return covariance(iBin,iBin-1);
0068     else return 0.;
0069   }
0070
0071   // TODO: verify this is actually correct
0072   double fullyCorrelatedUncertainty()
0073   {
0074     double sum = 0.;
0075
0076     for(int i=0;i<bins.size();i++)
0077     {
0078       for(int j=0;j<bins.size();j++)
0079       {
0080         if(i!=j) sum += covariance(i,j);
0081       }
0082     }
0083
0084     return sum/(pow(bins.size(),2)-bins.size());
0085   }
0086
0087   // TODO: implement
0088   double uncorrelatedUncertaintyInBin(double xlow, double xup)
0089   {
0090     std::cout << "WARNING: uncorrelatedUncertaintyInBin not yet implemented" << std::endl;
0091     return 0.;
0092   }
0093
0094   // TODO: implement
0095   double twoBinCovariance(double xlow, double xup, double xlow2, double xup2)
0096   {
0097     std::cout << "WARNING: twoBinCovariance not yet implemented" << std::endl;
0098     return 0.;
0099   }
0100
0101   protected:
0102
0103   TMatrixDSym covariance;
0104   std::vector<double> bins;
0105   std::unique_ptr<TH1D> hBins; // dummy histogram for easy bin-finding
0106
0107 };
0108
0109 #endif