diff --git a/inc/LauCalcChiSq.hh b/inc/LauCalcChiSq.hh --- a/inc/LauCalcChiSq.hh +++ b/inc/LauCalcChiSq.hh @@ -23,8 +23,12 @@ */ #include "TH2Poly.h" +#include "TGraph.h" #include "TString.h" +#include "ROOT/RDataFrame.hxx" +#include "ROOT/RResultPtr.hxx" + #include /*! \file LauCalcChiSq.hh @@ -115,6 +119,8 @@ TString treeName1_; //! Name of the high stats data tree TString treeName2_; + //! Experiment number in the toy file + TString iExpt_; //! Name of the x-coordinate branch in tree 1 TString xName1_; @@ -124,6 +130,10 @@ TString yName1_; //! Name of the y-coordinate branch in tree 2 TString yName2_; + //! Title of the x-coordinate for drawing + TString xTitle_; + //! Title of the y-coordinate for drawing + TString yTitle_; //! The minimum bin content Float_t minContent_; @@ -140,7 +150,13 @@ TH2Poly* chiSqHisto_; //! Histogram (constructed from template) filled with signed chisq of tree1 vs tree2 TH2Poly* chiSqSignedHisto_; - + //! Graph for data + ROOT::RDF::RResultPtr g_data_smartptr; + TGraph* g_data_; + //! Graph for toy + ROOT::RDF::RResultPtr g_toy_smartptr; + TGraph* g_toy_; + //! Minimum x coordinate of histograms Float_t xMin_; //! Maximum x coordinate of histograms diff --git a/src/LauCalcChiSq.cc b/src/LauCalcChiSq.cc --- a/src/LauCalcChiSq.cc +++ b/src/LauCalcChiSq.cc @@ -192,9 +192,14 @@ histo1_->SetTitleSize(0.06,"y"); histo1_->SetLabelSize(0.05,"x"); histo1_->SetLabelSize(0.05,"y"); - histo1_->SetXTitle(xName1_); - histo1_->SetYTitle(yName1_); - histo1_->Draw("colz"); + histo1_->SetMarkerColor(0); + g_data_->SetMarkerStyle(20); + 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"); histo2_->SetLabelFont(62,"x"); @@ -207,10 +212,20 @@ histo2_->SetLabelSize(0.05,"y"); histo2_->SetXTitle(xName1_); histo2_->SetYTitle(yName1_); - histo2_->Draw("colz"); - can.SaveAs("toy.pdf"); - - if(-1.*pullHisto_->GetMinimum() > pullHisto_->GetMaximum()) pullHisto_->SetMaximum(-1.*pullHisto_->GetMinimum()); + histo2_->SetMarkerColor(0); + g_toy_->SetMarkerStyle(20); + g_toy_->SetMarkerSize(.5); + 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()); pullHisto_->SetLabelFont(62,"x"); @@ -221,10 +236,20 @@ pullHisto_->SetTitleSize(0.06,"y"); pullHisto_->SetLabelSize(0.05,"x"); pullHisto_->SetLabelSize(0.05,"y"); - pullHisto_->SetXTitle(xName1_); - pullHisto_->SetYTitle(yName1_); - pullHisto_->Draw("colz"); + pullHisto_->SetMarkerColor(0); + pullHisto_->SetMinimum(-9.); + 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.C"); chiSqHisto_->SetLabelFont(62,"x"); chiSqHisto_->SetLabelFont(62,"y"); @@ -234,9 +259,10 @@ chiSqHisto_->SetTitleSize(0.06,"y"); chiSqHisto_->SetLabelSize(0.05,"x"); chiSqHisto_->SetLabelSize(0.05,"y"); - chiSqHisto_->SetXTitle(xName1_); - chiSqHisto_->SetYTitle(yName1_); - chiSqHisto_->Draw("colz"); + chiSqHisto_->SetMarkerColor(0); + g_data_->Draw("ap"); + //chiSqHisto_->Draw("colztextsame"); + chiSqHisto_->Draw("colzsame"); can.SaveAs("chiSq.pdf"); chiSqSignedHisto_->SetLabelFont(62,"x"); @@ -247,9 +273,12 @@ chiSqSignedHisto_->SetTitleSize(0.06,"y"); chiSqSignedHisto_->SetLabelSize(0.05,"x"); chiSqSignedHisto_->SetLabelSize(0.05,"y"); - chiSqSignedHisto_->SetXTitle(xName1_); - chiSqSignedHisto_->SetYTitle(yName1_); - chiSqSignedHisto_->Draw("colz"); + chiSqSignedHisto_->SetXTitle(xTitle_); + chiSqSignedHisto_->SetYTitle(yTitle_); + chiSqSignedHisto_->SetMarkerColor(0); + g_data_->Draw("ap"); + //chiSqSignedHisto_->Draw("colztextsame"); + chiSqSignedHisto_->Draw("colzsame"); can.SaveAs("chiSqSigned.pdf"); } @@ -257,32 +286,26 @@ void LauCalcChiSq::initialiseHistos() { - // Open the input control file: - // Low_stat_file_name Low_stat_tree_name 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 - // Min_bin_content N_free_params Low/high_stat_histo_ratio xMin xMax yMin yMax + // Open the input control file: + // 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 + // Min_bin_content N_free_params Low/high_stat_histo_ratio xMin xMax yMin yMax std::ifstream getData(inputFileName_.Data()); // 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()) { std::cerr<<"Error. Could not read first line of the input file "<Exit(EXIT_FAILURE); } - // open the file that contains the low stat histogram - TFile * file1 = TFile::Open(fileName1_.Data(), "read"); - - if (file1 == 0) {gSystem->Exit(EXIT_FAILURE);} - - // retrieve the low stat histogram - TTree* tree1 = dynamic_cast(file1->Get(treeName1_.Data())); - - if (tree1 == 0) { - std::cerr<<"Error. Could not find the tree "<Exit(EXIT_FAILURE); - } + // Open the low stat TTree (we'll be interested in one experiment only) + ROOT::RDataFrame d1(treeName1_.Data(), fileName1_.Data()); // get the info on the high stat histogram getData >> fileName2_ >> treeName2_ >> xName2_ >> yName2_; @@ -293,24 +316,13 @@ // if both histograms are in the same file then just retrieve the // high stat histogram from the file we already have open - TFile * file2(0); - TTree* tree2(0); + ROOT::RDataFrame d2=0; if ( fileName2_ == fileName1_ ) { - tree2 = dynamic_cast(file1->Get(treeName2_.Data())); + d2 = d1; } // otherwise open the other file and retrieve the high stat histogram from there else { - file2 = TFile::Open(fileName2_.Data(), "read"); - - if (file2 == 0) {gSystem->Exit(EXIT_FAILURE);} - - tree2 = dynamic_cast(file2->Get(treeName2_.Data())); - } - - if (tree2 == 0) { - std::cerr<<"Error. Could not find the tree "<Exit(EXIT_FAILURE); + d2 = ROOT::RDataFrame(treeName2_.Data(), fileName2_.Data()); } // get the info on the minimum content, number of parameters and scalefactor @@ -330,72 +342,50 @@ // close the text file getData.close(); - std::cout<<"Using the files and trees: "<SetBranchAddress(xName1_.Data(),&x); - tree1->SetBranchAddress(yName1_.Data(),&y); - - Int_t nEntries = tree1->GetEntries(); + std::cout <<"Using the files and trees: "<( g_data_smartptr.GetPtr() ); - Double_t* xs = new Double_t[nEntries]; - Double_t* ys = new Double_t[nEntries]; - - for ( Int_t i=0; i < nEntries; ++i ) { - tree1->GetEntry( i ); - xs[i] = x; - ys[i] = y; - } + // Get toy data + g_toy_smartptr = d2.Range(0,5000).Graph(xName1_.Data(),yName1_.Data()); + g_toy_ = g_toy_smartptr.GetPtr(); theHisto_ = new TH2Poly("theHisto_", "", xMin_, xMax_, yMin_, yMax_); //select the number of divisions to get us closest to minContent entries per bin std::vector 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_ - 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(theHisto_->Clone("histo1_")); histo2_ = dynamic_cast(theHisto_->Clone("histo2_")); pullHisto_ = dynamic_cast(theHisto_->Clone("pullHisto_")); chiSqHisto_ = dynamic_cast(theHisto_->Clone("chiSqHisto_")); chiSqSignedHisto_ = dynamic_cast(theHisto_->Clone("chiSqSignedHisto_")); - delete[] xs; - delete[] ys; - - //fill the two histograms from the trees - TString drawString1, drawString2, weightString2; - drawString1 += yName1_; - drawString1 += ":"; - drawString1 += xName1_; - drawString1 += ">>histo1_"; - drawString2 += yName2_; - drawString2 += ":"; - drawString2 += xName2_; - drawString2 += ">>histo2_"; - weightString2 += scaleFactor_; - - tree1->Draw(drawString1); - tree2->Draw(drawString2,weightString2); + + auto fillPoly = [this] ( const Double_t x, const Double_t y ) { + histo1_->Fill( x,y ); + }; + d1.Foreach( fillPoly, { xName1_.Data(), yName1_.Data()} ); + + ROOT::DisableImplicitMT(); // This is not thread-safe (thanks to TH2Poly) + auto fillPolyWeighted = [this] ( const Double_t x, const Double_t y ) { + histo2_->Fill( x,y,scaleFactor_ ); + }; + d2.Foreach( fillPolyWeighted, { xName1_.Data(), yName1_.Data()} ); histo1_->SetDirectory(0); histo2_->SetDirectory(0); pullHisto_->SetDirectory(0); chiSqHisto_->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& divisions)