CutUtils.h

//A collection of handy methods for doing cuts etc to make sure things are working

#include <TEventList.h>
#include <TH1.h>
#include <TROOT.h>

unsigned int errorCode = 0;
double Chi2Test(TH1* h1, TH1 *h2, Option_t *option, Int_t constraint);

//Save an eventlist based on cut and returns the number passing it
int count(TTree& t, string name, string cut) {
    string cmd = ">>" + name;
    t.Draw(cmd.c_str(), cut.c_str());
    TEventList *list = (TEventList*)gDirectory->Get(name.c_str()); 

    return list->GetN();
}

double cut(TTree& t, string name, string quantity, double meanCut, double rmsCut, double numberCut, bool prependAbs=true, string cut="1>0") {
    string histname = "h" + name;
    string cmd = quantity + ">>" + histname + "(100)";    
    t.Draw(cmd.c_str(), cut.c_str());
    TH1* hist = (TH1*)gDirectory->Get(histname.c_str());

    if (hist->GetMean() > meanCut) {
    ++errorCode;
    cout << "Failed mean cut for: " << quantity << " cut: " << meanCut << endl;
    }
    if (hist->GetRMS() > rmsCut) {
    ++errorCode;
    cout << "Failed rms cut for: " << quantity << " cut: " << meanCut << endl;
    }
    
    ostringstream oss;
    if (prependAbs) oss << "abs";
    oss << quantity << ">" << numberCut << " && " << cut;
    return count(t, name, oss.str())/(double)t.GetEntriesFast();
}

double compareShape(TTree& t, const string& name, const string& quantity1, const string& quantity2, const string& cut="1>0") {
    string histName1 = "h" + name + "1";
    string cmd1 = quantity1 + ">>" + histName1 + "(100)";
    t.Draw(cmd1.c_str(), cut.c_str());
    TH1* hist1 = (TH1*)gDirectory->Get(histName1.c_str());
    
    string histName2 = "h" + name + "2";
    string cmd2 = quantity2 + ">>" + histName2 + "(100)";
    t.Draw(cmd2.c_str(), cut.c_str());
    TH1* hist2 = (TH1*)gDirectory->Get(histName2.c_str());
    
    //Must be integer things!
    
    if (hist1->GetNbinsX() != hist2->GetNbinsX() || hist1->GetBinCenter(1)!=hist2->GetBinCenter(1) 
    || hist1->GetBinCenter(hist1->GetNbinsX())!=hist2->GetBinCenter(hist2->GetNbinsX())) {
    double lowedge = hist1->GetBinCenter(1)<hist2->GetBinCenter(1) ? hist1->GetBinCenter(1) : hist2->GetBinCenter(1);
    double highedge = hist1->GetBinCenter(hist1->GetNbinsX())>hist2->GetBinCenter(hist2->GetNbinsX())
              ? hist1->GetBinCenter(hist1->GetNbinsX()) : hist2->GetBinCenter(hist2->GetNbinsX());
    strstream oss1;
    oss1 << quantity1 << ">>" << histName1 << "(100," << lowedge << "," << highedge << ")";
    cmd1 = oss1.str();
    delete hist1;
    t.Draw(cmd1.c_str(), cut.c_str());
    hist1 = (TH1*)gDirectory->Get(histName1.c_str());
    
    strstream oss2;
    oss2 << quantity2 << ">>" << histName2 << "(100," << lowedge << "," << highedge << ")";
    cmd2 = oss2.str();
    delete hist2;
    t.Draw(cmd2.c_str(), cut.c_str());
    hist2 = (TH1*)gDirectory->Get(histName2.c_str());
    }
 
    if (hist1->GetEntries() == 0 || hist2->GetEntries()==0) return -1;
    
    //double kol = hist1->KolmogorovTest(hist2, "D");
    //cout << "Kolmogorov probability: " << (kol*100) << "%" << endl;
    return Chi2Test(hist1, hist2, "", 1);
}

double Chi2Test(TH1* h1, TH1 *h2, Option_t *option, Int_t constraint) {
  //The Chi2 (Pearson's) test for differences between h and this histogram. 
  //a small value of prob indicates a significant difference between the distributions
  //
  //if the data was collected in such a way that the number of entries
  //in the first histogram is necessarily equal to the number of entries
  //in the second, the parameter _constraint_ must be made 1. Default is 0.
  //any additional constraints on the data lower the number of degrees of freedom
  //(i.e. increase constraint to more positive values) in accordance with
  //their number
  //
  //"O" -overflows included
  //"U" - underflows included
  //
  //"P" - print information about number of degrees of freedom and
  //the value of chi2
  //by default underflows and overflows are not included

  //algorithm taken from "Numerical Recipes in C++"
  // implementation by Anna Kreshuk
  
  Int_t df;
  Double_t chsq = 0;
  Double_t prob;
  Double_t temp;
  Double_t koef1,koef2;
  Double_t nen1, nen2;
  Double_t bin1, bin2;
  Int_t i_start, i_end;

  TString opt = option;
  opt.ToUpper();

  TAxis *axis1 = h1->GetXaxis();
  TAxis *axis2 = h2->GetXaxis();

  Int_t nbins1 = axis1->GetNbins();
  Int_t nbins2 = axis2->GetNbins();

  //check dimensions
  if (h1->GetDimension()!=1 || h2->GetDimension()!=1){
     cout << "Chi2Test for 1-d only" << endl;
     return 0;
  }
  //check that the histograms are not empty
  nen1 = h1->GetEntries();
  nen2 = h2->GetEntries();
  if((nen1==0)||(nen2==0)){
     cout << "Chi2Test one of the histograms is empty" << endl;
     return 0;
  }
  //check number of channels
  if (nbins1 != nbins2){
     cout << "Chi2Test different number of channels" << endl;
     return 0;
  }

  //check binning
  Double_t diffprec = 1e-5;
  Double_t diff1 = TMath::Abs(axis1->GetXmin() - axis2->GetXmin());
  Double_t diff2 = TMath::Abs(axis1->GetXmax() - axis2->GetXmax());
  if ((diff1 > diffprec)||(diff2>diffprec)){
     cout << "Chi2Test different binning" << endl;
     return 0;
  }

  //see options

  i_start = 1;
  i_end = nbins1;
  df=nbins1-constraint;

  if(opt.Contains("O")) {
     i_end = nbins1+1;
     df++;
  }
  if(opt.Contains("U")) {
     i_start = 0;
     df++;
  }

   //the test
  if (TMath::Abs(nen1-nen2) > diffprec){
    koef1 = TMath::Sqrt(nen2/nen1);
    koef2 = TMath::Sqrt(nen1/nen2);
  } else{
    koef1 = 1;
    koef2 = 1;
  }

  for (Int_t i=i_start; i<=i_end; i++){
     bin1 = h1->GetBinContent(i);
     bin2 = h2->GetBinContent(i);
     if (bin1 ==0 && bin2==0){
        --df; //no data means one less degree of freedom
     } else {
        temp  = koef1*bin1-koef2*bin2;
        chsq += temp*temp/(bin1+bin2);
     }
  }
  
  prob = TMath::Prob(0.5*chsq, Int_t(0.5*df));
  
  if (opt.Contains("P")){
     cout << "the value of chi2 = " << chsq  << endl;
     cout << "the number of degrees of freedom = " << df << endl;
  }

  return prob;
}

---