72 inputFileName_(inputFileName),
102 std::cout<<
"Running chi-squared algorithm"<<std::endl;
105 std::cout<<
"Calculating chi-squared"<<std::endl;
112 std::cout<<
"Writing out histogram output"<<std::endl;
114 TFile* outFile =
new TFile(
"checkHistos.root",
"recreate");
116 histo1_->SetDirectory(outFile);
117 histo2_->SetDirectory(outFile);
125 std::cout<<
"Done"<<std::endl;
130 Float_t toyVal(0.), dataVal(0.);
131 Int_t minContent(0), ndof(0);
132 Float_t diff(0.), chiSq(0.), totalChiSq(0.);
134 minContent =
histo1_->GetBinContent(1);
136 for(Int_t i=1; i<=
histo1_->GetNumberOfBins(); ++i) {
138 if(
histo1_->GetBinContent(i)<minContent) {
139 minContent=
histo1_->GetBinContent(i);
145 toyVal =
histo2_->GetBinContent(i);
146 dataVal =
histo1_->GetBinContent(i);
147 diff = dataVal-toyVal;
148 if(toyVal>0) chiSq = (diff*diff)/toyVal;
164 std::cout<<
"Total ChiSq/nDof = "<<totalChiSq<<
"/"<<ndof<<
" = "<<totalChiSq/ndof<<std::endl;
165 std::cout <<
"Actual minimum entries per bin: " << minContent << std::endl;
170 gStyle->SetPalette(1,0);
171 const Int_t NRGBs = 5;
172 const Int_t NCont = 255;
173 Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00};
174 Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51};
175 Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00};
176 Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00};
177 TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
178 gStyle->SetNumberContours(NCont);
179 gStyle->SetOptStat(0000);
183 can.SetTopMargin(0.05);
184 can.SetRightMargin(0.17);
185 can.SetBottomMargin(0.16);
186 can.SetLeftMargin(0.14);
192 histo1_->SetTitleSize(0.06,
"x");
193 histo1_->SetTitleSize(0.06,
"y");
194 histo1_->SetLabelSize(0.05,
"x");
195 histo1_->SetLabelSize(0.05,
"y");
199 can.SaveAs(
"data.pdf");
205 histo2_->SetTitleSize(0.06,
"x");
206 histo2_->SetTitleSize(0.06,
"y");
207 histo2_->SetLabelSize(0.05,
"x");
208 histo2_->SetLabelSize(0.05,
"y");
212 can.SaveAs(
"toy.pdf");
228 can.SaveAs(
"pull.pdf");
241 can.SaveAs(
"chiSq.pdf");
254 can.SaveAs(
"chiSqSigned.pdf");
269 if (!getData.good()) {
270 std::cerr<<
"Error. Could not read first line of the input file "<<
inputFileName_<<std::endl;
271 gSystem->Exit(EXIT_FAILURE);
275 TFile * file1 = TFile::Open(
fileName1_.Data(),
"read");
277 if (file1 == 0) {gSystem->Exit(EXIT_FAILURE);}
280 TTree* tree1 =
dynamic_cast<TTree*
>(file1->Get(
treeName1_.Data()));
283 std::cerr<<
"Error. Could not find the tree "<<
treeName1_
285 gSystem->Exit(EXIT_FAILURE);
290 if (!getData.good()) {
291 std::cerr<<
"Error. Could not read the second line of the input file "<<
inputFileName_<<std::endl;
292 gSystem->Exit(EXIT_FAILURE);
300 tree2 =
dynamic_cast<TTree*
>(file1->Get(
treeName2_.Data()));
304 file2 = TFile::Open(
fileName2_.Data(),
"read");
306 if (file2 == 0) {gSystem->Exit(EXIT_FAILURE);}
308 tree2 =
dynamic_cast<TTree*
>(file2->Get(
treeName2_.Data()));
312 std::cerr<<
"Error. Could not find the tree "<<
treeName2_
314 gSystem->Exit(EXIT_FAILURE);
318 Int_t nParameters(0);
319 Float_t minContent(0.0), scaleFactor(1.0), xMin(0.0), xMax(1.0), yMin(0.0), yMax(1.0);
320 getData >> minContent >> nParameters >> scaleFactor >> xMin >> xMax >> yMin >> yMax;
321 if (getData.good()) {
336 std::cout<<
"Relative scaling factor histo1/histo2 set to "<<
scaleFactor_<<std::endl;
337 std::cout<<
"Minimum bin content is "<<
minContent_<<std::endl;
338 std::cout<<
"Number of free parameters is "<<
nParams_<<std::endl;
343 tree1->SetBranchAddress(
xName1_.Data(),&x);
344 tree1->SetBranchAddress(yName1_.Data(),&y);
346 Int_t nEntries = tree1->GetEntries();
348 Double_t* xs =
new Double_t[nEntries];
349 Double_t* ys =
new Double_t[nEntries];
351 for ( Int_t i=0; i < nEntries; ++i ) {
352 tree1->GetEntry( i );
360 std::vector<Int_t> divisions;
375 TString drawString1, drawString2, weightString2;
379 drawString1 +=
">>histo1_";
383 drawString2 +=
">>histo2_";
386 tree1->Draw(drawString1);
387 tree2->Draw(drawString2,weightString2);
396 if (file1 != 0) {file1->Close();}
398 if (file2 != 0) {file2->Close();}
406 for(Int_t i=0; i<nEntries; ++i) {
415 std::cout <<
"Target is " <<
minContent_ <<
" entries per bin" << std::endl;
416 std::cout <<
"Aiming to divide " << nIn <<
" entries between " << nBinsTarget <<
" bins" << std::endl;
420 Int_t nDivisions(0), nNines(0), nTwentyFives(0), nFortyNines(0), nEleventyElevens(0), nBins(1);
421 Int_t nDivisionsBest(0), nNinesBest(0), nTwentyFivesBest(0), nFortyNinesBest(0), nEleventyElevensBest(0), nBinsBest(1);
425 for(nNines=0; nNines<=nDivisions; ++nNines) {
426 for(nTwentyFives=0; nTwentyFives<=nDivisions-nNines; ++nTwentyFives) {
427 for(nFortyNines=0; nFortyNines<=nDivisions-nNines-nTwentyFives; ++nFortyNines) {
428 for(nEleventyElevens=0; nEleventyElevens<=nDivisions-nNines-nTwentyFives-nFortyNines; ++nEleventyElevens) {
429 nBins = TMath::Power(4,nDivisions-nNines-nTwentyFives-nFortyNines-nEleventyElevens)
430 *TMath::Power(9,nNines)*TMath::Power(25,nTwentyFives)
431 *TMath::Power(49,nFortyNines)*TMath::Power(121,nEleventyElevens);
432 if(nBins < nBinsTarget && nBins > nBinsBest) {
435 nDivisionsBest = nDivisions;
437 nTwentyFivesBest = nTwentyFives;
438 nFortyNinesBest = nFortyNines;
439 nEleventyElevensBest = nEleventyElevens;
445 }
while(TMath::Power(4,nDivisions+1) < nBinsTarget);
447 std::cout <<
"Using " << nBinsBest <<
" bins" << std::endl;
450 for(Int_t i=0; i<nEleventyElevensBest; ++i) {
451 divisions.push_back(11);
453 for(Int_t i=0; i<nFortyNinesBest; ++i) {
454 divisions.push_back(7);
456 for(Int_t i=0; i<nTwentyFivesBest; ++i) {
457 divisions.push_back(5);
459 for(Int_t i=0; i<nNinesBest; ++i) {
460 divisions.push_back(3);
462 for(Int_t i=0; i<nDivisionsBest-nNinesBest-nTwentyFivesBest-nFortyNinesBest-nEleventyElevensBest; ++i) {
463 divisions.push_back(2);
467 void LauCalcChiSq::getHisto(
const Double_t xMin,
const Double_t xMax,
const Double_t yMin,
const Double_t yMax,
const Double_t* xs,
const Double_t* ys,
const Int_t nEntries,
const std::vector<Int_t>& divisions,
const UInt_t iter)
471 if(iter == divisions.size()) {
472 Double_t * x_new =
new Double_t[5];
473 Double_t * y_new =
new Double_t[5];
474 x_new[0] = xMin; x_new[1] = xMin; x_new[2] = xMax; x_new[3] = xMax; x_new[4] = xMin;
475 y_new[0] = yMin; y_new[1] = yMax; y_new[2] = yMax; y_new[3] = yMin; y_new[4] = yMin;
477 if(
verbose_) std::cout <<
"Adding bin from (" << xMin <<
"," << yMin <<
") to (" << xMax <<
"," << yMax <<
")" << std::endl;
482 Int_t n_divx=divisions[iter];
483 Int_t n_divy=divisions[iter];
485 if(
verbose_) std::cout <<
"Dividing bin from (" << xMin <<
"," << yMin <<
") to (" << xMax <<
"," << yMax <<
") into " << n_divx <<
" by " << n_divy <<
" subbins" << std::endl;
487 Double_t *xIn =
new Double_t[nEntries];
488 Double_t *yIn =
new Double_t[nEntries];
489 Int_t *xIndex =
new Int_t [nEntries+2];
490 Int_t *yIndex =
new Int_t [nEntries+2];
493 for(Int_t i = 0; i<nEntries; ++i) {
494 if ((xs[i]<xMin)||(xs[i]>xMax)||(ys[i]<yMin)||(ys[i]>yMax))
continue;
495 xIn[xCountIn] = xs[i];
500 Double_t xLimits[n_divx + 1];
501 Double_t yLimits[n_divx][n_divy + 1];
504 TMath::Sort(xCountIn, xIn, xIndex,
false);
507 xLimits[n_divx] = xMax;
508 for (Int_t nDivx = 0; nDivx < n_divx; ++nDivx){
509 if (nDivx < (n_divx-1)){
510 xLimits[nDivx+1] = xIn[xIndex[xCountIn*(nDivx+1)/n_divx]];
514 yLimits[nDivx][0] = yMin;
515 yLimits[nDivx][n_divy] = yMax;
518 for(Int_t i = 0; i<nEntries; ++i) {
519 if ((xs[i]<xMin)||(xs[i]>xMax)||(ys[i]<yMin)||(ys[i]>yMax))
continue;
520 if ((xs[i]<xLimits[nDivx])||(xs[i]>=xLimits[nDivx+1])||(ys[i]<yMin)||(ys[i]>yMax))
continue;
521 yIn[yCountIn] = ys[i];
525 TMath::Sort(yCountIn, yIn, yIndex,
false);
527 for (Int_t nDivy = 1; nDivy < n_divy; ++nDivy){
528 yLimits[nDivx][nDivy] = yIn[yIndex[yCountIn*nDivy/n_divy]];
538 for (Int_t nDivx = 0; nDivx < n_divx; ++nDivx){
539 for (Int_t nDivy = 0; nDivy < n_divy; ++nDivy){
540 this->
getHisto(xLimits[nDivx], xLimits[nDivx + 1], yLimits[nDivx][nDivy], yLimits[nDivx][nDivy + 1], xs, ys, nEntries, divisions,iter+1);
TH2Poly * chiSqSignedHisto_
Histogram (constructed from template) filled with signed chisq of tree1 vs tree2. ...
TString treeName2_
Name of the high stats data tree.
TString xName1_
Name of the x-coordinate branch in tree 1.
File containing declaration of LauCalcChiSq class.
Bool_t verbose_
Verbose flag.
Int_t nParams_
Number of free parameters in fit (used for calculating the ndof)
void getHisto(const Double_t xMin, const Double_t xMax, const Double_t yMin, const Double_t yMax, const Double_t *xs, const Double_t *ys, const Int_t nEntries, const std::vector< Int_t > &divisions, const UInt_t iter=0)
Create the template histogram based on the binning scheme.
void calculateChiSq()
Calculate the chisq from the data histograms.
Float_t yMax_
Maximum y coordinate of histograms.
Float_t xMax_
Maximum x coordinate of histograms.
TString yName1_
Name of the y-coordinate branch in tree 1.
TString fileName2_
Name of the high stats data file.
TH2Poly * histo1_
Histogram (constructed from template) filled from tree 1.
TString fileName1_
Name of the low stats data file.
TH2Poly * pullHisto_
Histogram (constructed from template) filled with pulls of tree1 vs tree2.
void run()
Run the calculations.
Float_t yMin_
Minimum y coordinate of histograms.
Float_t xMin_
Minimum x coordinate of histograms.
TString xName2_
Name of the x-coordinate branch in tree 2.
TString inputFileName_
Name of the config file.
TH2Poly * theHisto_
Template histogram constructed from the binning scheme.
void pickBinning(const Double_t *xs, const Double_t *ys, const Int_t nEntries, std::vector< Int_t > &divisions)
Choose the binning scheme.
void initialiseHistos()
Read the config file, read the data and create histograms.
Utility class to allow the calculation of the chisq of the fit to the Dalitz plot.
TH2Poly * histo2_
Histogram (constructed from template) filled from tree 2.
TString treeName1_
Name of the low stats data tree.
Float_t scaleFactor_
Scalefactor between low and high stats data samples.
TH2Poly * chiSqHisto_
Histogram (constructed from template) filled with chisq of tree1 vs tree2.
virtual ~LauCalcChiSq()
Destructor.
TString yName2_
Name of the y-coordinate branch in tree 2.
Float_t minContent_
The minimum bin content.
void makePlots()
Create plots.