Changeset View
Changeset View
Standalone View
Standalone View
src/LauCalcChiSq.cc
Show First 20 Lines • Show All 186 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 low stat TTree (we'll be interested in one experiment only) | ||||
TFile * file1 = TFile::Open(fileName1_.Data(), "read"); | ROOT::RDataFrame d1(treeName1_.Data(), fileName1_.Data()); | ||||
if (file1 == 0) {gSystem->Exit(EXIT_FAILURE);} | |||||
// retrieve the low stat histogram | |||||
TTree* tree1 = dynamic_cast<TTree*>(file1->Get(treeName1_.Data())); | |||||
if (tree1 == 0) { | |||||
std::cerr<<"Error. Could not find the tree "<<treeName1_ | |||||
<<" in the file "<<fileName1_<<std::endl; | |||||
gSystem->Exit(EXIT_FAILURE); | |||||
} | |||||
// get the info on the high stat histogram | // get the info on the high stat histogram | ||||
getData >> fileName2_ >> treeName2_ >> xName2_ >> yName2_; | getData >> fileName2_ >> treeName2_ >> xName2_ >> yName2_; | ||||
if (!getData.good()) { | if (!getData.good()) { | ||||
std::cerr<<"Error. Could not read the second line of the input file "<<inputFileName_<<std::endl; | std::cerr<<"Error. Could not read the second line of the input file "<<inputFileName_<<std::endl; | ||||
gSystem->Exit(EXIT_FAILURE); | gSystem->Exit(EXIT_FAILURE); | ||||
} | } | ||||
// if both histograms are in the same file then just retrieve the | // if both histograms are in the same file then just retrieve the | ||||
// high stat histogram from the file we already have open | // high stat histogram from the file we already have open | ||||
TFile * file2(0); | ROOT::RDataFrame d2=0; | ||||
TTree* tree2(0); | |||||
if ( fileName2_ == fileName1_ ) { | if ( fileName2_ == fileName1_ ) { | ||||
tree2 = dynamic_cast<TTree*>(file1->Get(treeName2_.Data())); | d2 = d1; | ||||
} | } | ||||
// otherwise open the other file and retrieve the high stat histogram from there | // otherwise open the other file and retrieve the high stat histogram from there | ||||
else { | else { | ||||
file2 = TFile::Open(fileName2_.Data(), "read"); | d2 = ROOT::RDataFrame(treeName2_.Data(), fileName2_.Data()); | ||||
if (file2 == 0) {gSystem->Exit(EXIT_FAILURE);} | |||||
tree2 = dynamic_cast<TTree*>(file2->Get(treeName2_.Data())); | |||||
} | |||||
if (tree2 == 0) { | |||||
std::cerr<<"Error. Could not find the tree "<<treeName2_ | |||||
<<" in the file "<<fileName2_<<std::endl; | |||||
gSystem->Exit(EXIT_FAILURE); | |||||
} | } | ||||
// get the info on the minimum content, number of parameters and scalefactor | // get the info on the minimum content, number of parameters and scalefactor | ||||
Int_t nParameters(0); | Int_t nParameters(0); | ||||
Float_t minContent(0.0), scaleFactor(1.0), xMin(0.0), xMax(1.0), yMin(0.0), yMax(1.0); | Float_t minContent(0.0), scaleFactor(1.0), xMin(0.0), xMax(1.0), yMin(0.0), yMax(1.0); | ||||
getData >> minContent >> nParameters >> scaleFactor >> xMin >> xMax >> yMin >> yMax; | getData >> minContent >> nParameters >> scaleFactor >> xMin >> xMax >> yMin >> yMax; | ||||
if (getData.good()) { | if (getData.good()) { | ||||
minContent_ = minContent; | minContent_ = minContent; | ||||
nParams_ = nParameters; | nParams_ = nParameters; | ||||
scaleFactor_ = scaleFactor; | scaleFactor_ = scaleFactor; | ||||
xMin_ = xMin; | xMin_ = xMin; | ||||
xMax_ = xMax; | xMax_ = xMax; | ||||
yMin_ = yMin; | yMin_ = yMin; | ||||
yMax_ = yMax; | yMax_ = yMax; | ||||
} | } | ||||
// close the text file | // close the text file | ||||
getData.close(); | getData.close(); | ||||
std::cout<<"Using the files and trees: "<<fileName1_<<"; "<<treeName1_ | std::cout <<"Using the files and trees: "<<fileName1_<<"; "<<treeName1_ | ||||
<<" and "<<fileName2_<<"; "<<treeName2_<<std::endl; | <<" and "<<fileName2_<<"; "<<treeName2_<<std::endl; | ||||
std::cout<<"Relative scaling factor histo1/histo2 set to "<<scaleFactor_<<std::endl; | std::cout <<"Relative scaling factor histo1/histo2 set to "<<scaleFactor_<<std::endl; | ||||
std::cout<<"Minimum bin content is "<<minContent_<<std::endl; | std::cout <<"Minimum bin content is "<<minContent_<<std::endl; | ||||
std::cout<<"Number of free parameters is "<<nParams_<<std::endl; | std::cout <<"Number of free parameters is "<<nParams_<<std::endl; | ||||
Double_t x; | g_data_smartptr = d1.Filter(Form("iExpt==%s",iExpt_.Data())).Graph(xName1_.Data(),yName1_.Data()); | ||||
Double_t y; | g_data_ = dynamic_cast<TGraph*>( g_data_smartptr.GetPtr() ); | ||||
tree1->SetBranchAddress(xName1_.Data(),&x); | |||||
tree1->SetBranchAddress(yName1_.Data(),&y); | |||||
Int_t nEntries = tree1->GetEntries(); | |||||
Double_t* xs = new Double_t[nEntries]; | // Get toy data | ||||
Double_t* ys = new Double_t[nEntries]; | g_toy_smartptr = d2.Range(0,5000).Graph(xName1_.Data(),yName1_.Data()); | ||||
g_toy_ = g_toy_smartptr.GetPtr(); | |||||
for ( Int_t i=0; i < nEntries; ++i ) { | |||||
tree1->GetEntry( i ); | |||||
xs[i] = x; | |||||
ys[i] = y; | |||||
} | |||||
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(g_data_->GetX(),g_data_->GetY(),g_data_->GetN(),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_, g_data_->GetX(),g_data_->GetY(),g_data_->GetN(), divisions); | ||||
histo1_ = dynamic_cast<TH2Poly*>(theHisto_->Clone("histo1_")); | histo1_ = dynamic_cast<TH2Poly*>(theHisto_->Clone("histo1_")); | ||||
histo2_ = dynamic_cast<TH2Poly*>(theHisto_->Clone("histo2_")); | histo2_ = dynamic_cast<TH2Poly*>(theHisto_->Clone("histo2_")); | ||||
pullHisto_ = dynamic_cast<TH2Poly*>(theHisto_->Clone("pullHisto_")); | pullHisto_ = dynamic_cast<TH2Poly*>(theHisto_->Clone("pullHisto_")); | ||||
chiSqHisto_ = dynamic_cast<TH2Poly*>(theHisto_->Clone("chiSqHisto_")); | chiSqHisto_ = dynamic_cast<TH2Poly*>(theHisto_->Clone("chiSqHisto_")); | ||||
chiSqSignedHisto_ = dynamic_cast<TH2Poly*>(theHisto_->Clone("chiSqSignedHisto_")); | chiSqSignedHisto_ = dynamic_cast<TH2Poly*>(theHisto_->Clone("chiSqSignedHisto_")); | ||||
delete[] xs; | |||||
delete[] ys; | |||||
//fill the two histograms from the trees | auto fillPoly = [this] ( const Double_t x, const Double_t y ) { | ||||
TString drawString1, drawString2, weightString2; | histo1_->Fill( x,y ); | ||||
drawString1 += yName1_; | }; | ||||
drawString1 += ":"; | d1.Foreach( fillPoly, { xName1_.Data(), yName1_.Data()} ); | ||||
drawString1 += xName1_; | |||||
drawString1 += ">>histo1_"; | ROOT::DisableImplicitMT(); // This is not thread-safe (thanks to TH2Poly) | ||||
drawString2 += yName2_; | auto fillPolyWeighted = [this] ( const Double_t x, const Double_t y ) { | ||||
drawString2 += ":"; | histo2_->Fill( x,y,scaleFactor_ ); | ||||
drawString2 += xName2_; | }; | ||||
drawString2 += ">>histo2_"; | d2.Foreach( fillPolyWeighted, { xName1_.Data(), yName1_.Data()} ); | ||||
weightString2 += scaleFactor_; | |||||
tree1->Draw(drawString1); | |||||
tree2->Draw(drawString2,weightString2); | |||||
histo1_->SetDirectory(0); | histo1_->SetDirectory(0); | ||||
histo2_->SetDirectory(0); | histo2_->SetDirectory(0); | ||||
pullHisto_->SetDirectory(0); | pullHisto_->SetDirectory(0); | ||||
chiSqHisto_->SetDirectory(0); | chiSqHisto_->SetDirectory(0); | ||||
chiSqSignedHisto_->SetDirectory(0); | chiSqSignedHisto_->SetDirectory(0); | ||||
// close the file(s) containing the trees | |||||
if (file1 != 0) {file1->Close();} | |||||
delete file1; | |||||
if (file2 != 0) {file2->Close();} | |||||
delete file2; | |||||
} | } | ||||
void LauCalcChiSq::pickBinning(const Double_t* xs, const Double_t* ys, const Int_t nEntries, std::vector<Int_t>& divisions) | void LauCalcChiSq::pickBinning(const Double_t* xs, const Double_t* ys, const Int_t nEntries, std::vector<Int_t>& divisions) | ||||
{ | { | ||||
//first check how many events we have within the histogram limits | //first check how many events we have within the histogram limits | ||||
Int_t nIn(0); | Int_t nIn(0); | ||||
for(Int_t i=0; i<nEntries; ++i) { | for(Int_t i=0; i<nEntries; ++i) { | ||||
if(xs[i]<xMax_ && xs[i] >= xMin_ && ys[i]<yMax_ && ys[i] >= yMin_) { | if(xs[i]<xMax_ && xs[i] >= xMin_ && ys[i]<yMax_ && ys[i] >= yMin_) { | ||||
▲ Show 20 Lines • Show All 137 Lines • Show Last 20 Lines |