Changeset View
Changeset View
Standalone View
Standalone View
src/LauCalcChiSq.cc
Show All 16 Lines | |||||
/* | /* | ||||
Laura++ package authors: | Laura++ package authors: | ||||
John Back | John Back | ||||
Paul Harrison | Paul Harrison | ||||
Thomas Latham | Thomas Latham | ||||
*/ | */ | ||||
// Code to produce an adaptive binning scheme and calculate the 2D chi-square between | /*! \file LauCalcChiSq.cc | ||||
// two datasets (e.g. low-stat data and high-stat toyMC). | \brief File containing implementation of LauCalcChiSq class. | ||||
// To run the code, do "./CalcChiSq chiSqInput.txt", where chiSqInput.txt is | |||||
// an input control file, and contains the following lines: | Code to produce an adaptive binning scheme and calculate the 2D chi-square | ||||
// Low_stat_file_name Low_stat_histo_name | between two datasets (e.g. low-stat data and high-stat toyMC). | ||||
// High_stat_file_name High_stat_histo_name | |||||
// Min_bin_content N_free_params Low/high stat_histo_ratio | Note that the low and high stat histograms must have the same bin axes ranges | ||||
and number of bins. | |||||
// Note that the low and high stat histograms must have the same bin axes | |||||
// ranges and number of bins. | It works by using the low stat (first) histogram to find a binning scheme such | ||||
that the total number of entries in each bin is >= Min_bin_content. The number | |||||
// It works by using the low stat (first) histogram to find a binning scheme such | of entries in the histogram is divided by the desired minimum bin content to | ||||
// that the total number of entries in each bin is >= Min_bin_content. The number | give a target number of bins. The largest number of bins that can be expressed | ||||
// of entries in the histogram is divided by the desired minimum bin content to | as a product of powers of 4, 9, 25, 49 and 121 that does not exceed the target | ||||
// give a target number of bins. The largest number of bins that can be expressed | value is chosen. The histogram is the recursively subdivided in 2x2, 3x3, 5x5, | ||||
// as a product of powers of four, nine, 25, 49 and 121 that does not exceed the | 7x7 or 11x11 bins. For each stage of the subdivision, each bin is first | ||||
// target value is chosen. The histogram is the recursively subdivided in 2x2, 3x3, | divided into equally populated bins in x then each of these is further divded | ||||
// 5x5, 7x7 or 11x11 bins. For each stage of the subdivision, each bin is first | into equally popiulated bins in y. | ||||
// divided into equally populated bins in x then each of these is further divded | |||||
// into equally popiulated bins in y. | The (Pearson) chi-squared is then the sum of the chi-squared contributions of | ||||
all bins: | |||||
// The (Pearson) chi-squared is then the sum of the chi-squared contributions | (low_stat_number - high_stat_number)^2/(high_stat_number) | ||||
// of all bins: | The nDof = number of bins - number of free params - 1 | ||||
// (low_stat_number - high_stat_number)^2/(high_stat_number) | */ | ||||
// The nDof = number of bins - number of free params - 1 | |||||
#include "LauCalcChiSq.hh" | #include "LauCalcChiSq.hh" | ||||
#include "TAxis.h" | #include "TAxis.h" | ||||
#include "TFile.h" | #include "TFile.h" | ||||
#include "TMath.h" | #include "TMath.h" | ||||
#include "TSystem.h" | #include "TSystem.h" | ||||
#include "TTree.h" | #include "TTree.h" | ||||
▲ Show 20 Lines • Show All 130 Lines • ▼ Show 20 Lines | void LauCalcChiSq::makePlots() | ||||
histo1_->SetLabelFont(62,"x"); | histo1_->SetLabelFont(62,"x"); | ||||
histo1_->SetLabelFont(62,"y"); | histo1_->SetLabelFont(62,"y"); | ||||
histo1_->SetTitleFont(62,"x"); | histo1_->SetTitleFont(62,"x"); | ||||
histo1_->SetTitleFont(62,"y"); | histo1_->SetTitleFont(62,"y"); | ||||
histo1_->SetTitleSize(0.06,"x"); | histo1_->SetTitleSize(0.06,"x"); | ||||
histo1_->SetTitleSize(0.06,"y"); | histo1_->SetTitleSize(0.06,"y"); | ||||
histo1_->SetLabelSize(0.05,"x"); | histo1_->SetLabelSize(0.05,"x"); | ||||
histo1_->SetLabelSize(0.05,"y"); | histo1_->SetLabelSize(0.05,"y"); | ||||
histo1_->SetXTitle(xName1_); | histo1_->SetMarkerColor(0); | ||||
histo1_->SetYTitle(yName1_); | g_data_->SetMarkerStyle(20); | ||||
histo1_->Draw("colz"); | g_data_->SetMarkerSize(.5); | ||||
g_data_->Draw("ap"); | |||||
g_data_->GetXaxis()->SetTitle(xTitle_); | |||||
g_data_->GetYaxis()->SetTitle(yTitle_); | |||||
//histo1_->Draw("colztextsame"); | |||||
histo1_->Draw("colzsame"); | |||||
can.SaveAs("data.pdf"); | can.SaveAs("data.pdf"); | ||||
histo2_->SetLabelFont(62,"x"); | histo2_->SetLabelFont(62,"x"); | ||||
histo2_->SetLabelFont(62,"y"); | histo2_->SetLabelFont(62,"y"); | ||||
histo2_->SetTitleFont(62,"x"); | histo2_->SetTitleFont(62,"x"); | ||||
histo2_->SetTitleFont(62,"y"); | histo2_->SetTitleFont(62,"y"); | ||||
histo2_->SetTitleSize(0.06,"x"); | histo2_->SetTitleSize(0.06,"x"); | ||||
histo2_->SetTitleSize(0.06,"y"); | histo2_->SetTitleSize(0.06,"y"); | ||||
histo2_->SetLabelSize(0.05,"x"); | histo2_->SetLabelSize(0.05,"x"); | ||||
histo2_->SetLabelSize(0.05,"y"); | histo2_->SetLabelSize(0.05,"y"); | ||||
histo2_->SetXTitle(xName1_); | histo2_->SetXTitle(xName1_); | ||||
histo2_->SetYTitle(yName1_); | histo2_->SetYTitle(yName1_); | ||||
histo2_->Draw("colz"); | histo2_->SetMarkerColor(0); | ||||
can.SaveAs("toy.pdf"); | g_toy_->SetMarkerStyle(20); | ||||
g_toy_->SetMarkerSize(.5); | |||||
if(-1.*pullHisto_->GetMinimum() > pullHisto_->GetMaximum()) pullHisto_->SetMaximum(-1.*pullHisto_->GetMinimum()); | g_toy_->Draw("ap"); | ||||
g_toy_->GetXaxis()->SetTitle(xTitle_); | |||||
g_toy_->GetYaxis()->SetTitle(yTitle_); | |||||
//histo2_->Draw("colztextsame"); | |||||
histo2_->Draw("colzsame"); | |||||
can.SaveAs("toy.png"); | |||||
if(-1.*pullHisto_->GetMinimum() > pullHisto_->GetMaximum()) { | |||||
pullHisto_->SetMinimum(pullHisto_->GetMinimum()); // Required to workaround a bug in TH2Poly where, if SetMinimum is not called, negative entries are not drawn by "Draw("text")" | |||||
pullHisto_->SetMaximum(-1.*pullHisto_->GetMinimum()); | |||||
} | |||||
else pullHisto_->SetMinimum(-1.*pullHisto_->GetMaximum()); | else pullHisto_->SetMinimum(-1.*pullHisto_->GetMaximum()); | ||||
pullHisto_->SetLabelFont(62,"x"); | pullHisto_->SetLabelFont(62,"x"); | ||||
pullHisto_->SetLabelFont(62,"y"); | pullHisto_->SetLabelFont(62,"y"); | ||||
pullHisto_->SetTitleFont(62,"x"); | pullHisto_->SetTitleFont(62,"x"); | ||||
pullHisto_->SetTitleFont(62,"y"); | pullHisto_->SetTitleFont(62,"y"); | ||||
pullHisto_->SetTitleSize(0.06,"x"); | pullHisto_->SetTitleSize(0.06,"x"); | ||||
pullHisto_->SetTitleSize(0.06,"y"); | pullHisto_->SetTitleSize(0.06,"y"); | ||||
pullHisto_->SetLabelSize(0.05,"x"); | pullHisto_->SetLabelSize(0.05,"x"); | ||||
pullHisto_->SetLabelSize(0.05,"y"); | pullHisto_->SetLabelSize(0.05,"y"); | ||||
pullHisto_->SetXTitle(xName1_); | pullHisto_->SetMarkerColor(0); | ||||
pullHisto_->SetYTitle(yName1_); | pullHisto_->SetMinimum(-9.); | ||||
pullHisto_->Draw("colz"); | pullHisto_->SetMaximum(9.); | ||||
g_data_->SetTitle("(data - toy) / error"); | |||||
g_data_->Draw("ap"); | |||||
gStyle->SetPaintTextFormat(".2f"); | |||||
pullHisto_->SetLineWidth(1); | |||||
pullHisto_->SetLineColor(kBlack); | |||||
pullHisto_->SetLineStyle(1); | |||||
//pullHisto_->Draw("colztextsame"); | |||||
pullHisto_->Draw("colzsame"); | |||||
g_data_->Draw("psame"); | |||||
can.SaveAs("pull.pdf"); | can.SaveAs("pull.pdf"); | ||||
can.SaveAs("pull.C"); | |||||
chiSqHisto_->SetLabelFont(62,"x"); | chiSqHisto_->SetLabelFont(62,"x"); | ||||
chiSqHisto_->SetLabelFont(62,"y"); | chiSqHisto_->SetLabelFont(62,"y"); | ||||
chiSqHisto_->SetTitleFont(62,"x"); | chiSqHisto_->SetTitleFont(62,"x"); | ||||
chiSqHisto_->SetTitleFont(62,"y"); | chiSqHisto_->SetTitleFont(62,"y"); | ||||
chiSqHisto_->SetTitleSize(0.06,"x"); | chiSqHisto_->SetTitleSize(0.06,"x"); | ||||
chiSqHisto_->SetTitleSize(0.06,"y"); | chiSqHisto_->SetTitleSize(0.06,"y"); | ||||
chiSqHisto_->SetLabelSize(0.05,"x"); | chiSqHisto_->SetLabelSize(0.05,"x"); | ||||
chiSqHisto_->SetLabelSize(0.05,"y"); | chiSqHisto_->SetLabelSize(0.05,"y"); | ||||
chiSqHisto_->SetXTitle(xName1_); | chiSqHisto_->SetMarkerColor(0); | ||||
chiSqHisto_->SetYTitle(yName1_); | g_data_->Draw("ap"); | ||||
chiSqHisto_->Draw("colz"); | //chiSqHisto_->Draw("colztextsame"); | ||||
chiSqHisto_->Draw("colzsame"); | |||||
can.SaveAs("chiSq.pdf"); | can.SaveAs("chiSq.pdf"); | ||||
chiSqSignedHisto_->SetLabelFont(62,"x"); | chiSqSignedHisto_->SetLabelFont(62,"x"); | ||||
chiSqSignedHisto_->SetLabelFont(62,"y"); | chiSqSignedHisto_->SetLabelFont(62,"y"); | ||||
chiSqSignedHisto_->SetTitleFont(62,"x"); | chiSqSignedHisto_->SetTitleFont(62,"x"); | ||||
chiSqSignedHisto_->SetTitleFont(62,"y"); | chiSqSignedHisto_->SetTitleFont(62,"y"); | ||||
chiSqSignedHisto_->SetTitleSize(0.06,"x"); | chiSqSignedHisto_->SetTitleSize(0.06,"x"); | ||||
chiSqSignedHisto_->SetTitleSize(0.06,"y"); | chiSqSignedHisto_->SetTitleSize(0.06,"y"); | ||||
chiSqSignedHisto_->SetLabelSize(0.05,"x"); | chiSqSignedHisto_->SetLabelSize(0.05,"x"); | ||||
chiSqSignedHisto_->SetLabelSize(0.05,"y"); | chiSqSignedHisto_->SetLabelSize(0.05,"y"); | ||||
chiSqSignedHisto_->SetXTitle(xName1_); | chiSqSignedHisto_->SetXTitle(xTitle_); | ||||
chiSqSignedHisto_->SetYTitle(yName1_); | chiSqSignedHisto_->SetYTitle(yTitle_); | ||||
chiSqSignedHisto_->Draw("colz"); | chiSqSignedHisto_->SetMarkerColor(0); | ||||
g_data_->Draw("ap"); | |||||
//chiSqSignedHisto_->Draw("colztextsame"); | |||||
chiSqSignedHisto_->Draw("colzsame"); | |||||
can.SaveAs("chiSqSigned.pdf"); | can.SaveAs("chiSqSigned.pdf"); | ||||
} | } | ||||
void LauCalcChiSq::initialiseHistos() | void LauCalcChiSq::initialiseHistos() | ||||
{ | { | ||||
// Open the input control file: | // Open the input control file: | ||||
// Low_stat_file_name Low_stat_tree_name Low_stat_x_axis_name Low_stat_y_axis_name | // Low_stat_file_name Low_stat_tree_name iExpt Low_stat_x_axis_name Low_stat_y_axis_name | ||||
// High_stat_file_name High_stat_tree_name High_stat_x_axis_name High_stat_y_axis_name | // High_stat_file_name High_stat_tree_name High_stat_x_axis_name High_stat_y_axis_name | ||||
// Min_bin_content N_free_params Low/high_stat_histo_ratio xMin xMax yMin yMax | // Min_bin_content N_free_params Low/high_stat_histo_ratio xMin xMax yMin yMax | ||||
std::ifstream getData(inputFileName_.Data()); | std::ifstream getData(inputFileName_.Data()); | ||||
// get the info on the low stat histo | // get the info on the low stat histo | ||||
getData >> fileName1_ >> treeName1_ >> xName1_ >> yName1_; | getData >> fileName1_ >> treeName1_ >> iExpt_ >> xName1_ >> yName1_ >> xTitle_ >> yTitle_; | ||||
// Hack to introduce spaces in axis titles | |||||
xTitle_=xTitle_.ReplaceAll("SPACE"," "); | |||||
yTitle_=yTitle_.ReplaceAll("SPACE"," "); | |||||
if (!getData.good()) { | if (!getData.good()) { | ||||
std::cerr<<"Error. Could not read first line of the input file "<<inputFileName_<<std::endl; | std::cerr<<"Error. Could not read first line of the input file "<<inputFileName_<<std::endl; | ||||
gSystem->Exit(EXIT_FAILURE); | gSystem->Exit(EXIT_FAILURE); | ||||
} | } | ||||
// open the file that contains the low stat histogram | // open the file that contains the low stat histogram | ||||
TFile * file1 = TFile::Open(fileName1_.Data(), "read"); | TFile * file1 = TFile::Open(fileName1_.Data(), "read"); | ||||
if (file1 == 0) {gSystem->Exit(EXIT_FAILURE);} | if (file1 == 0) {gSystem->Exit(EXIT_FAILURE);} | ||||
// retrieve the low stat histogram | // retrieve the low stat histogram and extract the experiment | ||||
TTree* tree1 = dynamic_cast<TTree*>(file1->Get(treeName1_.Data())); | TTree* tree1temp = dynamic_cast<TTree*>(file1->Get(treeName1_.Data())); | ||||
TTree* tree1 = tree1temp->CopyTree("iExpt=="+iExpt_); | |||||
if (tree1 == 0) { | if (tree1 == 0) { | ||||
std::cerr<<"Error. Could not find the tree "<<treeName1_ | std::cerr<<"Error. Could not find the tree "<<treeName1_ | ||||
<<" in the file "<<fileName1_<<std::endl; | <<" in the file "<<fileName1_<<std::endl; | ||||
gSystem->Exit(EXIT_FAILURE); | gSystem->Exit(EXIT_FAILURE); | ||||
} | } | ||||
// get the info on the high stat histogram | // get the info on the high stat histogram | ||||
▲ Show 20 Lines • Show All 60 Lines • ▼ Show 20 Lines | void LauCalcChiSq::initialiseHistos() | ||||
Double_t* ys = new Double_t[nEntries]; | Double_t* ys = new Double_t[nEntries]; | ||||
for ( Int_t i=0; i < nEntries; ++i ) { | for ( Int_t i=0; i < nEntries; ++i ) { | ||||
tree1->GetEntry( i ); | tree1->GetEntry( i ); | ||||
xs[i] = x; | xs[i] = x; | ||||
ys[i] = y; | ys[i] = y; | ||||
} | } | ||||
g_data_ = new TGraph(nEntries,xs,ys); | |||||
// Get toy data | |||||
tree2->SetBranchAddress(xName1_.Data(),&x); | |||||
tree2->SetBranchAddress(yName1_.Data(),&y); | |||||
Int_t nEntries_toy = 5000;//tree2->GetEntries(); | |||||
Double_t* xs_toy = new Double_t[nEntries_toy]; | |||||
Double_t* ys_toy = new Double_t[nEntries_toy]; | |||||
for ( Int_t i=0; i < nEntries_toy; ++i ) { | |||||
tree2->GetEntry( i ); | |||||
xs_toy[i] = x; | |||||
ys_toy[i] = y; | |||||
} | |||||
g_toy_ = new TGraph(nEntries_toy,xs_toy,ys_toy); | |||||
theHisto_ = new TH2Poly("theHisto_", "", xMin_, xMax_, yMin_, yMax_); | theHisto_ = new TH2Poly("theHisto_", "", xMin_, xMax_, yMin_, yMax_); | ||||
//select the number of divisions to get us closest to minContent entries per bin | //select the number of divisions to get us closest to minContent entries per bin | ||||
std::vector<Int_t> divisions; | std::vector<Int_t> divisions; | ||||
this->pickBinning(xs,ys,nEntries,divisions); | this->pickBinning(xs,ys,nEntries,divisions); | ||||
//perform the adaptive bining based on histo1_ | //perform the adaptive bining based on histo1_ | ||||
this->getHisto(xMin_, xMax_, yMin_, yMax_, xs, ys, nEntries, divisions); | this->getHisto(xMin_, xMax_, yMin_, yMax_, xs, ys, nEntries, divisions); | ||||
▲ Show 20 Lines • Show All 180 Lines • Show Last 20 Lines |