LightFlavorRatios/corrections/CorrectionHistogram1D.h

0001 #ifndef CORRECTIONHISTOGRAM1D_H
0002 #define CORRECTIONHISTOGRAM1D_H
0003
0004 struct CorrectionHistogram1D
0005 {
0006   TH1* h_corr;
0007   std::string name;
0008   std::string title;
0009
0010   CorrectionHistogram1D() = default;
0011
0012   CorrectionHistogram1D(const CorrectionHistogram1D& ch)
0013   {
0014     h_corr = ch.h_corr;
0015     name = ch.name;
0016     title = ch.title;
0017   }
0018
0019   double get_xmin()
0020   {
0021     return h_corr->GetBinLowEdge(1);
0022   }
0023
0024   double get_xmax()
0025   {
0026     return h_corr->GetBinLowEdge(h_corr->GetNbinsX())+h_corr->GetBinWidth(h_corr->GetNbinsX());
0027   }
0028
0029   bool in_range(double x)
0030   {
0031     return (x>=get_xmin() && x<=get_xmax());
0032   }
0033
0034   void check_range_warn(double x_low, double x_high)
0035   {
0036     if(!in_range(x_low) || !in_range(x_high))
0037     {
0038       std::cout << "WARNING: selection [" << x_low << ", " << x_high << "]"
0039         << " is out of range [" << get_xmin() << ", "
0040         << get_xmax() << "]" << std::endl;
0041     }
0042   }
0043
0044   // get correction value by integrating appropriately over range of bins
0045   std::pair<double,double> get_val_and_error(double x_low, double x_high)
0046   {
0047     check_range_warn(x_low,x_high);
0048
0049     // add small buffer to avoid edge cases
0050     int bin_low = h_corr->FindBin(x_low+fabs(.0001*x_low));
0051     int bin_high = h_corr->FindBin(x_high-fabs(.0001*x_high));
0052
0053     if(bin_low == bin_high)
0054     {
0055       return std::make_pair(h_corr->GetBinContent(bin_low),h_corr->GetBinError(bin_low));
0056     }
0057     else
0058     {
0059       double binwidth_low = h_corr->GetBinWidth(bin_low);
0060       double binwidth_high = h_corr->GetBinWidth(bin_high);
0061
0062       double upperedge_low = h_corr->GetBinLowEdge(bin_low)+binwidth_low;
0063       double loweredge_high = h_corr->GetBinLowEdge(bin_high);
0064
0065       double endbin_frac_low = (upperedge_low - x_low)/binwidth_low;
0066       double endbin_frac_high = (x_high - loweredge_high)/binwidth_high;
0067
0068       // weighted average by bin error
0069       double weighted_sum = 0.;
0070       double sum_inv_w2 = 0.;
0071
0072       // add ends separately
0073       weighted_sum += h_corr->GetBinContent(bin_low)*endbin_frac_low/pow(h_corr->GetBinError(bin_low),2.);
0074       sum_inv_w2 += endbin_frac_low/pow(h_corr->GetBinError(bin_low),2.);
0075
0076       weighted_sum += h_corr->GetBinContent(bin_high)*endbin_frac_high/pow(h_corr->GetBinError(bin_high),2.);
0077       sum_inv_w2 += endbin_frac_high/pow(h_corr->GetBinError(bin_high),2.);
0078
0079       // and everything else in between
0080       for(int ibin = bin_low+1; ibin<bin_high; ibin++)
0081       {
0082         weighted_sum += h_corr->GetBinContent(ibin)/pow(h_corr->GetBinError(ibin),2.);
0083         sum_inv_w2 += 1./pow(h_corr->GetBinError(ibin),2.);
0084       }
0085
0086       return std::make_pair(weighted_sum/sum_inv_w2,1/sqrt(sum_inv_w2));
0087     }
0088   }
0089
0090   virtual void apply_correction(float xlow, float xhigh, TH1F* h, int bin)
0091   {
0092     std::cout << "WARNING: you haven't replaced apply_correction with a correction-specific function yet!!" << std::endl;
0093   }
0094
0095   virtual ~CorrectionHistogram1D() = default;
0096 };
0097
0098 #endif