LightFlavorRatios/util/HistogramTools.h

0001 #ifndef HISTOGRAM_TOOLS_H
0002 #define HISTOGRAM_TOOLS_H
0003
0004 #include <TH1F.h>
0005 #include <TH2F.h>
0006 #include <TTree.h>
0007
0008 #include <regex>
0009
0010 std::vector<std::string> cutvars_from_cutstring(std::string cutstring)
0011 {
0012   // all possible tokens dividing cut variables
0013   std::regex dividers("[&|!<=>]+");
0014   // strip whitespace
0015   cutstring.erase(std::remove_if(cutstring.begin(),cutstring.end(),::isspace),cutstring.end());
0016   std::vector<std::string> statement_list(std::sregex_token_iterator(cutstring.begin(),cutstring.end(),dividers,-1),
0017                                           std::sregex_token_iterator());
0018   std::vector<std::string> cutvars;
0019   for(std::string statement : statement_list)
0020   {
0021     // is this a numeric value?
0022     bool is_numeric;
0023     try
0024     {
0025       size_t nparsed = 0;
0026       std::stod(statement,&nparsed);
0027       is_numeric = (nparsed == statement.size());
0028     }
0029     catch(const std::invalid_argument&)
0030     {
0031       is_numeric = false;
0032     }
0033     catch(const std::out_of_range&)
0034     {
0035       is_numeric = false;
0036     }
0037     // cut variables are the non-numeric strings
0038     if(!is_numeric)
0039     {
0040       cutvars.push_back(statement);
0041     }
0042   }
0043
0044   return cutvars;
0045 }
0046
0047 struct HistogramInfo
0048 {
0049   std::string name;
0050   std::string title;
0051   std::string axis_label;
0052   std::vector<double> bins;
0053
0054   HistogramInfo(const std::string& hname, const std::string& htitle, const std::vector<double>& hbins,
0055                 const std::string& haxislabel = "")
0056   : name(hname), title(htitle), axis_label(haxislabel), bins(hbins)
0057   {
0058   }
0059
0060   // constructor for uniform binning
0061   HistogramInfo(const std::string& hname, const std::string& htitle, const int nBins, const float xmin, const float xmax,
0062                 const std::string& haxislabel = "")
0063   : name(hname), title(htitle), axis_label(haxislabel)
0064   {
0065     const float interval = (xmax-xmin)/nBins;
0066     for(int i=0; i<=nBins; i++)
0067     {
0068       bins.push_back(xmin+i*interval);
0069     }
0070   }
0071
0072   std::string get_bin_selection(std::string var, size_t bin) const
0073   {
0074     if(bin==0) // underflow
0075     {
0076       return var+"<"+std::to_string(bins[0]);
0077     }
0078     else if(bin<bins.size()) // Nbins = bins.size()-1
0079     {
0080       return var+">="+std::to_string(bins[bin-1])+"&&"+var+"<"+std::to_string(bins[bin]);
0081     }
0082     else if(bin==bins.size()) // overflow
0083     {
0084       return var+">="+std::to_string(bins[bins.size()-1]);
0085     }
0086     else
0087     {
0088       return "";
0089     }
0090   }
0091 };
0092
0093 std::vector<std::string> get_cutvars(TTree* t,std::string cut_string)
0094 {
0095   std::vector<std::string> cut_vars;
0096
0097   TObjArray* branches = t->GetListOfBranches();
0098   size_t nbranches = t->GetNbranches();
0099   for(size_t i=0;i<nbranches;i++)
0100   {
0101     std::string branchname = branches->At(i)->GetName();
0102     if(cut_string.find(branchname) != std::string::npos)
0103     {
0104       cut_vars.push_back(branchname);
0105     }
0106   }
0107
0108   return cut_vars;
0109 }
0110
0111 bool isIntBranch(TBranch* b)
0112 {
0113   TClass* c;
0114   EDataType type;
0115   b->GetExpectedType(c,type);
0116   std::vector<EDataType> integral_datatypes = {
0117     kInt_t,
0118     kUInt_t,
0119     kChar_t,
0120     kUChar_t,
0121     kShort_t,
0122     kUShort_t,
0123     kLong_t,
0124     kULong_t,
0125     kLong64_t,
0126     kULong64_t
0127   };
0128   return std::find(integral_datatypes.begin(),integral_datatypes.end(),type)!=integral_datatypes.end();
0129 }
0130
0131 int findBin(float val, std::vector<double> bins)
0132 {
0133   int thisbin = -1;
0134   for(size_t i=0; i<bins.size(); i++)
0135   {
0136     if(val<bins[i])
0137     {
0138       thisbin = i;
0139       break;
0140     }
0141   }
0142   // if val is not < any bin edge, it belongs in last bin
0143   if(thisbin == -1)
0144   {
0145     thisbin = bins.size();
0146   }
0147   return thisbin;
0148 }
0149
0150 std::vector<double> makeLogBins(const int nBins, const float xmin, const float xmax)
0151 {
0152   std::vector<double> bins;
0153   const float log_interval = (log(xmax)-log(xmin))/nBins;
0154   for(int i=0; i<=nBins; i++)
0155   {
0156     bins.push_back(exp(log(xmin) + i*log_interval));
0157   }
0158   return bins;
0159 }
0160
0161 TH1F* makeHistogram(const HistogramInfo& hinfo)
0162 {
0163   TH1F* h = new TH1F(hinfo.name.c_str(),(hinfo.title+";"+hinfo.axis_label+";Candidates").c_str(),hinfo.bins.size()-1,hinfo.bins.data());
0164   return h;
0165 }
0166
0167 TH1F* makeHistogram(const std::string& basename, const std::string& basetitle, const HistogramInfo& hinfo)
0168 {
0169   TH1F* h = new TH1F((basename+"_vs"+hinfo.name).c_str(),(basetitle+" vs. "+hinfo.title+";"+hinfo.axis_label+";Candidates").c_str(),hinfo.bins.size()-1,hinfo.bins.data());
0170   return h;
0171 }
0172
0173 TH2F* make2DHistogram(const std::string& basename, const std::string& basetitle, const HistogramInfo& hinfo_x, const HistogramInfo& hinfo_y)
0174 {
0175   const std::string name = basename + hinfo_x.name + "_vs" + hinfo_y.name;
0176   const std::string title = basetitle + " " + hinfo_x.title + " vs. " + hinfo_y.title + ";" + hinfo_x.axis_label + ";" + hinfo_y.axis_label;
0177   TH2F* h = new TH2F(name.c_str(),title.c_str(),hinfo_x.bins.size()-1,hinfo_x.bins.data(),hinfo_y.bins.size()-1,hinfo_y.bins.data());
0178   return h;
0179 }
0180
0181 std::vector<TH1F*> makeDifferentialHistograms(const HistogramInfo& hinfo_x, const HistogramInfo& hinfo_y)
0182 {
0183   std::vector<TH1F*> h_out;
0184   for(size_t i=0; i<=hinfo_y.bins.size(); i++)
0185   {
0186     const std::string name = hinfo_x.name + "_vs" + hinfo_y.name + "_" + std::to_string(i);
0187     const std::string title = hinfo_x.title + " vs. " + hinfo_y.title +
0188       " (bin " + std::to_string(i) + ", " + hinfo_y.get_bin_selection(hinfo_y.name,i) + ");" + hinfo_x.axis_label;
0189     h_out.push_back(new TH1F(name.c_str(),title.c_str(),hinfo_x.bins.size()-1,hinfo_x.bins.data()));
0190   }
0191   return h_out;
0192 }
0193
0194 #endif // HISTOGRAM_TOOLS_H