diff --git a/src/LauCalcChiSq.cc b/src/LauCalcChiSq.cc
index 5d7e31c..b687c25 100644
--- a/src/LauCalcChiSq.cc
+++ b/src/LauCalcChiSq.cc
@@ -1,544 +1,543 @@
 
 /*
 Copyright 2008 University of Warwick
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at
 
     http://www.apache.org/licenses/LICENSE-2.0
 
 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.
 */
 
 /*
 Laura++ package authors:
 John Back
 Paul Harrison
 Thomas Latham
 */
 
-// Code to produce an adaptive binning scheme and calculate the 2D chi-square between
-// two datasets (e.g. low-stat data and high-stat toyMC).
-// To run the code, do "./CalcChiSq chiSqInput.txt", where chiSqInput.txt is
-// an input control file, and contains the following lines:
-// Low_stat_file_name      Low_stat_histo_name
-// 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.
-
-// 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
-// of entries in the histogram is divided by the desired minimum bin content to
-// give a target number of bins. The largest number of bins that can be expressed
-// as a product of powers of four, nine, 25, 49 and 121 that does not exceed the
-// target value is chosen. The histogram is the recursively subdivided in 2x2, 3x3,
-// 5x5, 7x7 or 11x11 bins. For each stage of the subdivision, each bin is first
-// 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:
-// (low_stat_number - high_stat_number)^2/(high_stat_number)
-// The nDof = number of bins - number of free params - 1
+/*! \file LauCalcChiSq.cc
+  \brief File containing implementation of LauCalcChiSq class.
+
+  Code to produce an adaptive binning scheme and calculate the 2D chi-square
+  between two datasets (e.g. low-stat data and high-stat toyMC).
+ 
+  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
+  of entries in the histogram is divided by the desired minimum bin content to
+  give a target number of bins. The largest number of bins that can be expressed
+  as a product of powers of 4, 9, 25, 49 and 121 that does not exceed the target
+  value is chosen. The histogram is the recursively subdivided in 2x2, 3x3, 5x5,
+  7x7 or 11x11 bins. For each stage of the subdivision, each bin is first
+  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:
+  (low_stat_number - high_stat_number)^2/(high_stat_number)
+  The nDof = number of bins - number of free params - 1
+ */
 
 #include "LauCalcChiSq.hh"
 
 #include "TAxis.h"
 #include "TFile.h"
 #include "TMath.h"
 #include "TSystem.h"
 #include "TTree.h"
 
 #include "TCanvas.h"
 #include "TColor.h"
 #include "TStyle.h"
 
 #include <iostream>
 #include <fstream>
 #include <cstdlib>
 #include <cmath>
 
 ClassImp(LauCalcChiSq)
 
 
 LauCalcChiSq::LauCalcChiSq(const TString& inputFileName) :
 	inputFileName_(inputFileName),
 	fileName1_(""),
 	fileName2_(""),
 	treeName1_(""),
 	treeName2_(""),
 	xName1_(""),
 	xName2_(""),
 	yName1_(""),
 	yName2_(""),
 	minContent_(10.0),
 	histo1_(0),
 	histo2_(0),
 	chiSqHisto_(0),
 	chiSqSignedHisto_(0),
 	xMin_(0.0),
 	xMax_(0.0),
 	yMin_(0.0),
 	yMax_(0.0),
 	nParams_(0),
 	scaleFactor_(1.0),
 	verbose_(kFALSE)
 {
 }
 
 LauCalcChiSq::~LauCalcChiSq()
 {
 }
 
 void LauCalcChiSq::run()
 {
 	std::cout<<"Running chi-squared algorithm"<<std::endl;
 	this->initialiseHistos();
 
 	std::cout<<"Calculating chi-squared"<<std::endl;
 	this->calculateChiSq();
 
 	//make plots
 	this->makePlots();
 
 	// Output the various histograms
 	std::cout<<"Writing out histogram output"<<std::endl;
 
 	TFile* outFile = new TFile("checkHistos.root", "recreate");
 
 	histo1_->SetDirectory(outFile);
 	histo2_->SetDirectory(outFile);
 	pullHisto_->SetDirectory(outFile);
 	chiSqHisto_->SetDirectory(outFile);
 	chiSqSignedHisto_->SetDirectory(outFile);
 
 	outFile->Write();
 	outFile->Close();
 
 	std::cout<<"Done"<<std::endl;
 }
 
 void LauCalcChiSq::calculateChiSq()
 {
 	Float_t toyVal(0.), dataVal(0.);
 	Int_t minContent(0), ndof(0);
 	Float_t diff(0.), chiSq(0.), totalChiSq(0.);
 
 	minContent = histo1_->GetBinContent(1);
 
 	for(Int_t i=1; i<=histo1_->GetNumberOfBins(); ++i) {
 		//keep track of actual minimum
 		if(histo1_->GetBinContent(i)<minContent) {
 			minContent=histo1_->GetBinContent(i);
 		}
 
 		// Calculate Pearson chi-square for this bin, using the
 		// second histogram for the expected distribution
 		chiSq = 0.;
 		toyVal  = histo2_->GetBinContent(i);
 		dataVal = histo1_->GetBinContent(i);
 		diff = dataVal-toyVal;
 		if(toyVal>0) chiSq = (diff*diff)/toyVal;
 		totalChiSq += chiSq;
 
 		chiSqHisto_->SetBinContent(i, chiSq);
 
 		if(diff>0) {
 			chiSqSignedHisto_->SetBinContent(i, chiSq);
 			pullHisto_->SetBinContent(i, sqrt(chiSq));
 		} else {
 			chiSqSignedHisto_->SetBinContent(i, -chiSq);
 			pullHisto_->SetBinContent(i, -sqrt(chiSq));
 		}
 	}
 
 	ndof = histo1_->GetNumberOfBins()-nParams_-1;
 
 	std::cout<<"Total ChiSq/nDof = "<<totalChiSq<<"/"<<ndof<<" = "<<totalChiSq/ndof<<std::endl;
 	std::cout << "Actual minimum entries per bin: " << minContent << std::endl;
 }
 
 void LauCalcChiSq::makePlots()
 {
    gStyle->SetPalette(1,0);
    const Int_t NRGBs = 5;
    const Int_t NCont = 255;
    Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00};
    Double_t red[NRGBs]   = { 0.00, 0.00, 0.87, 1.00, 0.51};
    Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00};
    Double_t blue[NRGBs]  = { 0.51, 1.00, 0.12, 0.00, 0.00};
    TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
    gStyle->SetNumberContours(NCont);
    gStyle->SetOptStat(0000);
 
    TCanvas can;
 
    can.SetTopMargin(0.05);
    can.SetRightMargin(0.17);
    can.SetBottomMargin(0.16);
    can.SetLeftMargin(0.14);
 
    histo1_->SetLabelFont(62,"x");
    histo1_->SetLabelFont(62,"y");
    histo1_->SetTitleFont(62,"x");
    histo1_->SetTitleFont(62,"y");
    histo1_->SetTitleSize(0.06,"x");
    histo1_->SetTitleSize(0.06,"y");
    histo1_->SetLabelSize(0.05,"x");
    histo1_->SetLabelSize(0.05,"y");
    histo1_->SetXTitle(xName1_);
    histo1_->SetYTitle(yName1_);
    histo1_->Draw("colz");
    can.SaveAs("data.pdf");
 
    histo2_->SetLabelFont(62,"x");
    histo2_->SetLabelFont(62,"y");
    histo2_->SetTitleFont(62,"x");
    histo2_->SetTitleFont(62,"y");
    histo2_->SetTitleSize(0.06,"x");
    histo2_->SetTitleSize(0.06,"y");
    histo2_->SetLabelSize(0.05,"x");
    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());
    else pullHisto_->SetMinimum(-1.*pullHisto_->GetMaximum());
 
    pullHisto_->SetLabelFont(62,"x");
    pullHisto_->SetLabelFont(62,"y");
    pullHisto_->SetTitleFont(62,"x");
    pullHisto_->SetTitleFont(62,"y");
    pullHisto_->SetTitleSize(0.06,"x");
    pullHisto_->SetTitleSize(0.06,"y");
    pullHisto_->SetLabelSize(0.05,"x");
    pullHisto_->SetLabelSize(0.05,"y");
    pullHisto_->SetXTitle(xName1_);
    pullHisto_->SetYTitle(yName1_);
    pullHisto_->Draw("colz");
    can.SaveAs("pull.pdf");
 
    chiSqHisto_->SetLabelFont(62,"x");
    chiSqHisto_->SetLabelFont(62,"y");
    chiSqHisto_->SetTitleFont(62,"x");
    chiSqHisto_->SetTitleFont(62,"y");
    chiSqHisto_->SetTitleSize(0.06,"x");
    chiSqHisto_->SetTitleSize(0.06,"y");
    chiSqHisto_->SetLabelSize(0.05,"x");
    chiSqHisto_->SetLabelSize(0.05,"y");
    chiSqHisto_->SetXTitle(xName1_);
    chiSqHisto_->SetYTitle(yName1_);
    chiSqHisto_->Draw("colz");
    can.SaveAs("chiSq.pdf");
 
    chiSqSignedHisto_->SetLabelFont(62,"x");
    chiSqSignedHisto_->SetLabelFont(62,"y");
    chiSqSignedHisto_->SetTitleFont(62,"x");
    chiSqSignedHisto_->SetTitleFont(62,"y");
    chiSqSignedHisto_->SetTitleSize(0.06,"x");
    chiSqSignedHisto_->SetTitleSize(0.06,"y");
    chiSqSignedHisto_->SetLabelSize(0.05,"x");
    chiSqSignedHisto_->SetLabelSize(0.05,"y");
    chiSqSignedHisto_->SetXTitle(xName1_);
    chiSqSignedHisto_->SetYTitle(yName1_);
    chiSqSignedHisto_->Draw("colz");
    can.SaveAs("chiSqSigned.pdf");
 
 }
 
 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
 	std::ifstream getData(inputFileName_.Data());
 
 	// get the info on the low stat histo
 	getData >> fileName1_ >> treeName1_ >> xName1_ >> yName1_;
 	if (!getData.good()) {
 		std::cerr<<"Error. Could not read first line of the input file "<<inputFileName_<<std::endl;
 		gSystem->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<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
 	getData >> fileName2_ >> treeName2_ >> xName2_ >> yName2_;
 	if (!getData.good()) {
 		std::cerr<<"Error. Could not read the second line of the input file "<<inputFileName_<<std::endl;
 		gSystem->Exit(EXIT_FAILURE);
 	}
 
 	// 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);
 	if ( fileName2_ == fileName1_ ) {
 		tree2 = dynamic_cast<TTree*>(file1->Get(treeName2_.Data()));
 	}
 	// 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<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
 	Int_t nParameters(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;
 	if (getData.good()) {
 		minContent_ = minContent;
 		nParams_ = nParameters;
 		scaleFactor_ = scaleFactor;
 		xMin_ = xMin;
 		xMax_ = xMax;
 		yMin_ = yMin;
 		yMax_ = yMax;
 	}
 
 	// close the text file
 	getData.close();
 
 	std::cout<<"Using the files and trees: "<<fileName1_<<"; "<<treeName1_
 	    <<" and "<<fileName2_<<"; "<<treeName2_<<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<<"Number of free parameters is "<<nParams_<<std::endl;
 
 	Double_t x;
 	Double_t y;
 
 	tree1->SetBranchAddress(xName1_.Data(),&x);
 	tree1->SetBranchAddress(yName1_.Data(),&y);
 
 	Int_t nEntries = tree1->GetEntries();
 
 	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;
 	}
 
 	theHisto_ = new TH2Poly("theHisto_", "", xMin_, xMax_, yMin_, yMax_);
 
 	//select the number of divisions to get us closest to minContent entries per bin
 	std::vector<Int_t> divisions;
 	this->pickBinning(xs,ys,nEntries,divisions);
 
 	//perform the adaptive bining based on  histo1_
 	this->getHisto(xMin_, xMax_, yMin_, yMax_, xs, ys, nEntries, divisions);
 
 	histo1_           = dynamic_cast<TH2Poly*>(theHisto_->Clone("histo1_"));
 	histo2_           = dynamic_cast<TH2Poly*>(theHisto_->Clone("histo2_"));
 	pullHisto_        = dynamic_cast<TH2Poly*>(theHisto_->Clone("pullHisto_"));
 	chiSqHisto_       = dynamic_cast<TH2Poly*>(theHisto_->Clone("chiSqHisto_"));
 	chiSqSignedHisto_ = dynamic_cast<TH2Poly*>(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);
 
 	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<Int_t>& divisions)
 {
 	//first check how many events we have within the histogram limits
 	Int_t nIn(0);
 	for(Int_t i=0; i<nEntries; ++i) {
 		if(xs[i]<xMax_ && xs[i] >= xMin_ && ys[i]<yMax_ && ys[i] >= yMin_) {
 			++nIn;
 		}
 	}
 
 	//aim to have exactly minContent events in each bin
 	Int_t nBinsTarget = nIn / minContent_;
 
 	std::cout << "Target is " << minContent_ << " entries per bin" << std::endl;
 	std::cout << "Aiming to divide " << nIn << " entries between " << nBinsTarget << " bins" << std::endl;
 
 	//we will iteratively sub-divide histogram bins into either 4, 9, 25, 49 or 121
 	//here we figure out how many 4s, 9s, 25s, 49s and 121s to use to best match our target without exceeding it
 	Int_t nDivisions(0), nNines(0), nTwentyFives(0), nFortyNines(0), nEleventyElevens(0), nBins(1);
 	Int_t nDivisionsBest(0), nNinesBest(0), nTwentyFivesBest(0), nFortyNinesBest(0), nEleventyElevensBest(0), nBinsBest(1);
 
 	do {
 		++nDivisions;
 		for(nNines=0; nNines<=nDivisions; ++nNines) {
 			for(nTwentyFives=0; nTwentyFives<=nDivisions-nNines; ++nTwentyFives) {
 				for(nFortyNines=0; nFortyNines<=nDivisions-nNines-nTwentyFives; ++nFortyNines) {
 					for(nEleventyElevens=0; nEleventyElevens<=nDivisions-nNines-nTwentyFives-nFortyNines; ++nEleventyElevens) {
 						nBins =  TMath::Power(4,nDivisions-nNines-nTwentyFives-nFortyNines-nEleventyElevens)
 							*TMath::Power(9,nNines)*TMath::Power(25,nTwentyFives)
 							*TMath::Power(49,nFortyNines)*TMath::Power(121,nEleventyElevens);
 						if(nBins < nBinsTarget && nBins > nBinsBest) {
 							//keep track of the best number of bins so far
 							nBinsBest = nBins;
 							nDivisionsBest = nDivisions;
 							nNinesBest = nNines;
 							nTwentyFivesBest = nTwentyFives;
 							nFortyNinesBest = nFortyNines;
 							nEleventyElevensBest = nEleventyElevens;
 						}
 					}
 				}
 			}
 		}
 	} while(TMath::Power(4,nDivisions+1) < nBinsTarget);//if 4^n > target then we've gone far enough
 
 	std::cout << "Using " << nBinsBest << " bins" << std::endl;
 
 	//fill the vector with the divisions that we want to make
 	for(Int_t i=0; i<nEleventyElevensBest; ++i) {
 		divisions.push_back(11);
 	}
 	for(Int_t i=0; i<nFortyNinesBest; ++i) {
 		divisions.push_back(7);
 	}
 	for(Int_t i=0; i<nTwentyFivesBest; ++i) {
 		divisions.push_back(5);
 	}
 	for(Int_t i=0; i<nNinesBest; ++i) {
 		divisions.push_back(3);
 	}
 	for(Int_t i=0; i<nDivisionsBest-nNinesBest-nTwentyFivesBest-nFortyNinesBest-nEleventyElevensBest; ++i) {
 		divisions.push_back(2);
 	}
 }
 
 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)
 {
 
 	//If it's the last iteration create the bin and return
 	if(iter == divisions.size()) {
 		Double_t * x_new = new Double_t[5];
 		Double_t * y_new = new Double_t[5];
 		x_new[0] = xMin; x_new[1] = xMin; x_new[2] = xMax; x_new[3] = xMax; x_new[4] = xMin;
 		y_new[0] = yMin; y_new[1] = yMax; y_new[2] = yMax; y_new[3] = yMin; y_new[4] = yMin;
 		theHisto_->AddBin(5, x_new, y_new);
 		if(verbose_) std::cout << "Adding bin from (" << xMin << "," << yMin << ") to (" << xMax << "," << yMax << ")" << std::endl;
 		return;
 	}
 
 	//If not the last iteration then divide the bin
 	Int_t n_divx=divisions[iter];
 	Int_t n_divy=divisions[iter];
 
 	if(verbose_) std::cout << "Dividing bin from (" << xMin << "," << yMin << ") to (" << xMax << "," << yMax << ") into " << n_divx << " by " << n_divy << " subbins" << std::endl;
 
 	Double_t *xIn = new Double_t[nEntries];
 	Double_t *yIn = new Double_t[nEntries];
 	Int_t *xIndex = new Int_t [nEntries+2];
 	Int_t *yIndex = new Int_t [nEntries+2];
 
 	Int_t xCountIn = 0;
 	for(Int_t i = 0; i<nEntries; ++i) {
 		if ((xs[i]<xMin)||(xs[i]>xMax)||(ys[i]<yMin)||(ys[i]>yMax)) continue;
 		xIn[xCountIn] = xs[i];
 		++xCountIn;
 	}
 
 	//find the delimitting x and y values for the sub bins
 	Double_t xLimits[n_divx + 1];
 	Double_t yLimits[n_divx][n_divy + 1];
 
 	//first sort entries in x and divide bin into equally populated bins in x
 	TMath::Sort(xCountIn, xIn, xIndex, false);
 
 	xLimits[0] = xMin;
 	xLimits[n_divx] = xMax;
 	for (Int_t nDivx = 0; nDivx < n_divx; ++nDivx){
 		if (nDivx  < (n_divx-1)){
 			xLimits[nDivx+1] = xIn[xIndex[xCountIn*(nDivx+1)/n_divx]];
 		}
 
 		//for each bin in x divide into equally populated bins in y
 		yLimits[nDivx][0] = yMin;
 		yLimits[nDivx][n_divy] = yMax;
 		Int_t yCountIn = 0;
 
 		for(Int_t i = 0; i<nEntries; ++i) {
 			if ((xs[i]<xMin)||(xs[i]>xMax)||(ys[i]<yMin)||(ys[i]>yMax)) continue;
 			if ((xs[i]<xLimits[nDivx])||(xs[i]>=xLimits[nDivx+1])||(ys[i]<yMin)||(ys[i]>yMax)) continue;
 			yIn[yCountIn] = ys[i];
 			++yCountIn;
 		}
 
 		TMath::Sort(yCountIn, yIn, yIndex, false);
 
 		for (Int_t nDivy = 1; nDivy < n_divy; ++nDivy){
 			yLimits[nDivx][nDivy] = yIn[yIndex[yCountIn*nDivy/n_divy]];
 		}
 	}
 
 	delete[] xIn;
 	delete[] yIn;
 	delete[] xIndex;
 	delete[] yIndex;
 
 	//call for each sub bin
 	for (Int_t nDivx = 0; nDivx < n_divx; ++nDivx){
 		for (Int_t nDivy = 0; nDivy < n_divy; ++nDivy){
 			this->getHisto(xLimits[nDivx], xLimits[nDivx + 1], yLimits[nDivx][nDivy], yLimits[nDivx][nDivy + 1], xs, ys, nEntries, divisions,iter+1);
 		}
 	}
 }
 
diff --git a/src/LauDabbaRes.cc b/src/LauDabbaRes.cc
index 764a05c..48e2b55 100644
--- a/src/LauDabbaRes.cc
+++ b/src/LauDabbaRes.cc
@@ -1,267 +1,268 @@
 
 /*
 Copyright 2010 University of Warwick
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at
 
     http://www.apache.org/licenses/LICENSE-2.0
 
 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.
 */
 
 /*
 Laura++ package authors:
 John Back
 Paul Harrison
 Thomas Latham
 */
 
 /*! \file LauDabbaRes.cc
     \brief File containing implementation of LauDabbaRes class.
+
     Formulae and data values from arXiv:0901.2217 - author D.V.Bugg
 */
 
 #include <iostream>
 
 #include "LauConstants.hh"
 #include "LauDabbaRes.hh"
 #include "LauResonanceInfo.hh"
 
 ClassImp(LauDabbaRes)
 
 
 LauDabbaRes::LauDabbaRes(LauResonanceInfo* resInfo, const Int_t resPairAmpInt, const LauDaughters* daughters) :
 	LauAbsResonance(resInfo, resPairAmpInt, daughters),
 	mSumSq_(0.0),
 	sAdler_(0.0),
 	b_(0),
 	alpha_(0),
 	beta_(0)
 {
 	// Default constant factors
 	const Double_t bVal     = 24.49;
 	const Double_t alphaVal = 0.1;
 	const Double_t betaVal  = 0.1;
 
 	const TString& parNameBase = this->getSanitisedName();
 
 	TString bName(parNameBase);
 	bName += "_b";
 	b_ = resInfo->getExtraParameter( bName );
 	if ( b_ == 0 ) {
 		b_ = new LauParameter( bName, bVal, 0.0, 100.0, kTRUE );
 		b_->secondStage(kTRUE);
 		resInfo->addExtraParameter( b_ );
 	}
 
 	TString alphaName(parNameBase);
 	alphaName += "_alpha";
 	alpha_ = resInfo->getExtraParameter( alphaName );
 	if ( alpha_ == 0 ) {
 		alpha_ = new LauParameter( alphaName, alphaVal, 0.0, 10.0, kTRUE );
 		alpha_->secondStage(kTRUE);
 		resInfo->addExtraParameter( alpha_ );
 	}
 
 	TString betaName(parNameBase);
 	betaName += "_beta";
 	beta_ = resInfo->getExtraParameter( betaName );
 	if ( beta_ == 0 ) {
 		beta_ = new LauParameter( betaName, betaVal, 0.0, 10.0, kTRUE );
 		beta_->secondStage(kTRUE);
 		resInfo->addExtraParameter( beta_ );
 	}
 }
 
 LauDabbaRes::~LauDabbaRes()
 {
 }
 
 void LauDabbaRes::initialise()
 {
 	// check that we have a D and a pi
 	this->checkDaughterTypes();
 
 	// Initialise various constants
 	Double_t massDaug1 = this->getMassDaug1();
 	Double_t massDaug2 = this->getMassDaug2();
 
 	Double_t mSum = massDaug1 + massDaug2;
 	mSumSq_ = mSum*mSum;
 
 	Double_t massDaug1Sq = massDaug1*massDaug1;
 	Double_t massDaug2Sq = massDaug2*massDaug2;
 	sAdler_ = TMath::Max(massDaug1Sq,massDaug2Sq) - 0.5*TMath::Min(massDaug1Sq,massDaug2Sq); // Adler zero at (mD)^2 - 0.5*(mpi)^2
 
 	Int_t resSpin = this->getSpin();
 	if (resSpin != 0) {
 		std::cerr << "WARNING in LauDabbaRes::initialise : Spin = " << resSpin << " is not zero!  It will be ignored anyway!" << std::endl;
 	}
 }
 
 void LauDabbaRes::checkDaughterTypes() const
 {
 	// Check that the daughter tracks are D and pi. Otherwise issue a warning.
 	Int_t resPairAmpInt = this->getPairInt();
 	if (resPairAmpInt < 1 || resPairAmpInt > 3) {
 		std::cerr << "WARNING in LauDabbaRes::checkDaughterTypes : resPairAmpInt = " << resPairAmpInt << " is out of the range [1,2,3]." << std::endl;
 		return;
 	}
 
 	// Check that daughter types agree
 	const TString& nameDaug1 = this->getNameDaug1();
 	const TString& nameDaug2 = this->getNameDaug2();
 	if ( !( nameDaug1.Contains("pi", TString::kIgnoreCase) && nameDaug2.Contains("d", TString::kIgnoreCase) ) ) {
 		if ( !( nameDaug2.Contains("pi", TString::kIgnoreCase) && nameDaug1.Contains("d", TString::kIgnoreCase) ) ) {
 			std::cerr << "ERROR in LauDabbaRes::checkDaughterTypes : Dabba model is using daughters \"" << nameDaug1 << "\" and \"" << nameDaug2 << "\" that are not a D meson and a pion." << std::endl;
 		}
 	}
 }
 
 LauComplex LauDabbaRes::resAmp(Double_t mass, Double_t spinTerm)
 {
 	// This function returns the complex dynamical amplitude for a Dabba distribution
 	// given the invariant mass and cos(helicity) values.
 
 	// Invariant mass squared combination for the system
 	Double_t s = mass*mass;
 
 	// Dabba is spin zero - so there are no helicity factors.
 	// Just set it to 1.0 in case anyone decides to use it at a later date.
 	spinTerm = 1.0;
 
 	// Phase-space factor
 	Double_t rho(0.0);
 	Double_t sDiff = s - mSumSq_;
 	if ( sDiff > 0.0 ) {
 		rho = TMath::Sqrt(1.0 - mSumSq_/s);
 	}
 
 	const Double_t bVal = this->getBValue();
 	const Double_t alphaVal = this->getAlphaValue();
 	const Double_t betaVal = this->getBetaValue();
 
 	Double_t realPart = 1.0 - betaVal * sDiff;
 
 	Double_t imagPart = bVal * TMath::Exp( - alphaVal * sDiff ) * ( s - sAdler_ ) * rho;
 
 	LauComplex resAmplitude( realPart, imagPart );
 
 	Double_t denomFactor = realPart*realPart + imagPart*imagPart;
 
 	Double_t invDenomFactor = 0.0;
 	if (denomFactor > 1e-10) {invDenomFactor = 1.0/denomFactor;}
 
 	resAmplitude.rescale(spinTerm*invDenomFactor);
 
 	return resAmplitude;
 }
 
 const std::vector<LauParameter*>& LauDabbaRes::getFloatingParameters()
 {
 	this->clearFloatingParameters();
 
 	if ( ! this->fixBValue() ) {
 		this->addFloatingParameter( b_ );
 	}
 
 	if ( ! this->fixAlphaValue() ) {
 		this->addFloatingParameter( alpha_ );
 	}
 
 	if ( ! this->fixBetaValue() ) {
 		this->addFloatingParameter( beta_ );
 	}
 
 	return this->getParameters();
 }
 
 void LauDabbaRes::setResonanceParameter(const TString& name, const Double_t value) 
 {
 	// Set various parameters for the lineshape
 	if (name == "b") {
 		this->setBValue(value);
 		std::cout << "INFO in LauDabbaRes::setResonanceParameter : Setting parameter b = " << this->getBValue() << std::endl;
 	}
 	else if (name == "alpha") {
 		this->setAlphaValue(value);
 		std::cout << "INFO in LauDabbaRes::setResonanceParameter : Setting parameter alpha = " << this->getAlphaValue() << std::endl;
 	}
 	else if (name == "beta") {
 		this->setBetaValue(value);
 		std::cout << "INFO in LauDabbaRes::setResonanceParameter : Setting parameter beta = " << this->getBetaValue() << std::endl;
 	}
 	else {
 		std::cerr << "WARNING in LauDabbaRes::setResonanceParameter: Parameter name not reconised.  No parameter changes made." << std::endl;
 	}
 }
 
 void LauDabbaRes::floatResonanceParameter(const TString& name)
 {
 	if (name == "b") {
 		if ( b_->fixed() ) {
 			b_->fixed( kFALSE );
 			this->addFloatingParameter( b_ );
 		} else {
 			std::cerr << "WARNING in LauDabbaRes::floatResonanceParameter: Parameter already floating.  No parameter changes made." << std::endl;
 		}
 	} else if (name == "alpha") {
 		if ( alpha_->fixed() ) {
 			alpha_->fixed( kFALSE );
 			this->addFloatingParameter( alpha_ );
 		} else {
 			std::cerr << "WARNING in LauDabbaRes::floatResonanceParameter: Parameter already floating.  No parameter changes made." << std::endl;
 		}
 	} else if (name == "beta") {
 		if ( beta_->fixed() ) {
 			beta_->fixed( kFALSE );
 			this->addFloatingParameter( beta_ );
 		} else {
 			std::cerr << "WARNING in LauDabbaRes::floatResonanceParameter: Parameter already floating.  No parameter changes made." << std::endl;
 		}
 	} else {
 		std::cerr << "WARNING in LauDabbaRes::fixResonanceParameter: Parameter name not reconised.  No parameter changes made." << std::endl;
 	}
 }
 
 LauParameter* LauDabbaRes::getResonanceParameter(const TString& name)
 {
 	if (name == "b") {
 		return b_;
 	} else if (name == "alpha") {
 		return alpha_;
 	} else if (name == "beta") {
 		return beta_;
 	} else {
 		std::cerr << "WARNING in LauDabbaRes::getResonanceParameter: Parameter name not reconised." << std::endl;
 		return 0;
 	}
 }
 
 void LauDabbaRes::setBValue(const Double_t b)
 {
 	b_->value( b );
 	b_->genValue( b );
 	b_->initValue( b );
 }
 
 void LauDabbaRes::setAlphaValue(const Double_t alpha)
 {
 	alpha_->value( alpha );
 	alpha_->genValue( alpha );
 	alpha_->initValue( alpha );
 }
 
 void LauDabbaRes::setBetaValue(const Double_t beta)
 {
 	beta_->value( beta );
 	beta_->genValue( beta );
 	beta_->initValue( beta );
 }
 
diff --git a/src/LauMergeDataFiles.cc b/src/LauMergeDataFiles.cc
index 46fa1ae..10328b1 100644
--- a/src/LauMergeDataFiles.cc
+++ b/src/LauMergeDataFiles.cc
@@ -1,306 +1,310 @@
 
 /*
 Copyright 2008 University of Warwick
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at
 
     http://www.apache.org/licenses/LICENSE-2.0
 
 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.
 */
 
 /*
 Laura++ package authors:
 John Back
 Paul Harrison
 Thomas Latham
 */
 
+/*! \file LauMergeDataFiles.cc
+  \brief File containing implementation of LauMergeDataFiles class.
+ */
+
 #include "LauMergeDataFiles.hh"
 
 #include <iostream>
 #include <map>
 
 #include "TLeaf.h"
 #include "TObjArray.h"
 #include "TSystem.h"
 
 ClassImp(LauMergeDataFiles)
 
 
 LauMergeDataFiles::LauMergeDataFiles(const TString& fileName1, const TString& fileName2, const TString& treeName) :
 	fileName1_(fileName1),
 	fileName2_(fileName2),
 	treeName_(treeName),
 	inputFile1_(0),
 	inputFile2_(0),
 	inputTree1_(0),
 	inputTree2_(0),
 	outputFile_(0),
 	outputTree_(0)
 {
 }
 
 LauMergeDataFiles::~LauMergeDataFiles()
 {
 	if (inputFile1_ && inputFile1_->IsOpen()) {
 		inputFile1_->Close();
 	}
 	delete inputFile1_;
 
 	if (inputFile2_ && inputFile2_->IsOpen()) {
 		inputFile2_->Close();
 	}
 	delete inputFile2_;
 
 	if (outputFile_ && outputFile_->IsOpen()) {
 		outputFile_->Close();
 	}
 	delete outputFile_;
 }
 
 void LauMergeDataFiles::openInputFiles()
 {
 	// open the two input ROOT files
 	inputFile1_ = TFile::Open(fileName1_);
 	if (!inputFile1_ || inputFile1_->IsZombie()) {
 		std::cerr<<"Problem opening file: \""<<fileName1_<<"\", exiting..."<<std::endl;
 		gSystem->Exit(EXIT_FAILURE);
 	}
 	inputTree1_ = dynamic_cast<TTree*>( inputFile1_->Get(treeName_) );
 	if (!inputTree1_) {
 		std::cerr<<"Problem getting tree called "<<treeName_<<"from file: \""<<fileName1_<<"\", exiting..."<<std::endl;
 		gSystem->Exit(EXIT_FAILURE);
 	}
 
 	inputFile2_ = TFile::Open(fileName2_);
 	if (!inputFile2_ || inputFile2_->IsZombie()) {
 		std::cerr<<"Problem opening file: \""<<fileName2_<<"\", exiting..."<<std::endl;
 		gSystem->Exit(EXIT_FAILURE);
 	}
 	inputTree2_ = dynamic_cast<TTree*>( inputFile2_->Get(treeName_) );
 	if (!inputTree2_) {
 		std::cerr<<"Problem getting tree called "<<treeName_<<"from file: \""<<fileName2_<<"\", exiting..."<<std::endl;
 		gSystem->Exit(EXIT_FAILURE);
 	}
 }
 
 void LauMergeDataFiles::setupInputTrees()
 {
 	TObjArray* leaves1 = inputTree1_->GetListOfLeaves();
 	TObjArray* leaves2 = inputTree2_->GetListOfLeaves();
 	Int_t nLeaves1 = leaves1->GetEntries();
 	Int_t nLeaves2 = leaves2->GetEntries();
 	if ( nLeaves1 != nLeaves2 ) {
 		std::cerr<<"Different number of leaves in the two input trees, not continuing."<<std::endl;
 		return;
 	}
 
 	std::cout<<"Setting branches for input trees \""<<treeName_<<"\"..."<<std::endl;
 	inputTree1_->SetBranchAddress("iExpt",&iExpt_);
 	inputTree1_->SetBranchAddress("iEvtWithinExpt",&iEvtWithinExpt_);
 	inputTree2_->SetBranchAddress("iExpt",&iExpt_);
 	inputTree2_->SetBranchAddress("iEvtWithinExpt",&iEvtWithinExpt_);
 
 	for (Int_t iLeaf(0); iLeaf<nLeaves1; ++iLeaf) {
 
 		TLeaf * leaf = dynamic_cast<TLeaf*>((*leaves1)[iLeaf]);
 		TString type = leaf->GetTypeName();
 		TString name = leaf->GetName();
 		Int_t size = leaf->GetNdata();
 
 		if ((name == "iExpt") || (name == "iEvtWithinExpt") || (size != 1)) {
 			continue;
 		}
 
 		if ( type == "Double_t" ) {
 			std::pair<LeafDoubleMap::iterator,bool> result = doubleVars_.insert(std::make_pair(name,0.0));
 			LeafDoubleMap::iterator iter = result.first;
 			bool ok = result.second;
 			if (ok) {
 				inputTree1_->SetBranchAddress(name,&(iter->second));
 				inputTree2_->SetBranchAddress(name,&(iter->second));
 			}
 		} else if ( type == "Int_t" ) {
 			std::pair<LeafIntegerMap::iterator,bool> result = integerVars_.insert(std::make_pair(name,0));
 			LeafIntegerMap::iterator iter = result.first;
 			bool ok = result.second;
 			if (ok) {
 				inputTree1_->SetBranchAddress(name,&(iter->second));
 				inputTree2_->SetBranchAddress(name,&(iter->second));
 			}
 		}
 	}
 	std::cout<<"Set branch addresses for "<<doubleVars_.size()<<" Double_t branches.\n";
 	std::cout<<"Set branch addresses for "<<integerVars_.size()<<" Int_t branches.\n"<<std::endl;
 }
 
 void LauMergeDataFiles::setupOutputTree()
 {
 	std::cout<<"Creating branches for output tree \""<<outputTree_->GetName()<<"\"..."<<std::endl;
 
 	outputTree_->Branch("iExpt",&iExpt_,"iExpt/I");
 	outputTree_->Branch("iEvtWithinExpt",&iEvtWithinExpt_,"iEvtWithinExpt/I");
 
 	for (LeafDoubleMap::iterator iter = doubleVars_.begin(); iter != doubleVars_.end(); ++iter) {
 		TString name = iter->first;
 		Double_t * address = &(iter->second);
 		TString thirdBit = name;
 		thirdBit += "/D";
 
 		outputTree_->Branch(name,address,thirdBit);
 	}
 	for (LeafIntegerMap::iterator iter = integerVars_.begin(); iter != integerVars_.end(); ++iter) {
 		TString name = iter->first;
 		Int_t * address = &(iter->second);
 		TString thirdBit = name;
 		thirdBit += "/I";
 
 		outputTree_->Branch(name,address,thirdBit);
 	}
 	std::cout<<"Created "<<doubleVars_.size()+integerVars_.size()<<" branches.\n"<<std::endl;
 }
 
 void LauMergeDataFiles::process(const TString& fileName)
 {
 	this->openInputFiles();
 
 	this->setupInputTrees();
 
 	outputFile_ = TFile::Open(fileName,"recreate");
 	outputTree_ = new TTree(treeName_,"");
 	this->setupOutputTree();
 
 	// loop over the trees and combine the corresponding experiments
 	std::cout<<"Starting to combine the trees..."<<std::endl;
 
 	// Find the first and last entries for each experiment in each tree
 	this->findExperiments( inputTree1_, tree1Expts_ );
 	this->findExperiments( inputTree2_, tree2Expts_ );
 
 	// Check that the experiments in the two trees match
 	if ( !this->checkExperimentMaps() ) {
 		return;
 	}
 
 	// Loop through the experiments
 	for ( ExptsMap::const_iterator iter1 = tree1Expts_.begin(); iter1 != tree1Expts_.end(); ++iter1 ) {
 
 		// get the map element for tree2
 		Int_t expt = iter1->first;
 		ExptsMap::const_iterator iter2 = tree2Expts_.find( expt );
 
 		// determine the number of entries in tree1
 		Int_t nEntriesInTree1 = iter1->second.second - iter1->second.first + 1;
 
 		// read the entries from the trees, filling the output tree
 		this->readExperiment( inputTree1_, iter1, 0 );
 		this->readExperiment( inputTree2_, iter2, nEntriesInTree1 );
 	}
 
 	// Write the output file
 	this->writeFile();
 }
 
 void LauMergeDataFiles::findExperiments(TTree* tree, ExptsMap& exptsMap)
 {
 	const Int_t nEntries = tree->GetEntries();
 
 	// loop through the tree
 	for ( Int_t iEntry(0); iEntry<nEntries; ++iEntry ) {
 		// grab the entry
 		tree->GetEntry(iEntry);
 
 		// see if we already have an element in the map for the
 		// current experiment
 		ExptsMap::iterator iter = exptsMap.find(iExpt_);
 		if ( iter == exptsMap.end() ) {
 			// if not, we need to add an element that points to
 			// this entry in the tree as the start entry
 			exptsMap.insert( std::make_pair( iExpt_, std::make_pair( iEntry, -99 ) ) );
 
 			// also we need to complete the map element for the
 			// previous experiment with the previous tree entry
 			// as the last entry
 			ExptsMap::iterator previter = exptsMap.find(iExpt_-1);
 			if ( previter != exptsMap.end() ) {
 				previter->second.second = iEntry-1;
 			}
 		}
 	}
 	// need to complete the map element for the final experiment
 	exptsMap[iExpt_].second = nEntries-1;
 }
 
 Bool_t LauMergeDataFiles::checkExperimentMaps() const
 {
 	// first check that the two maps are the same size
 	UInt_t size1 = tree1Expts_.size();
 	UInt_t size2 = tree2Expts_.size();
 	if ( size1 != size2 ) {
 		std::cerr<<"ERROR in LauMergeDataFiles::checkExperimentMaps : Experiment maps are not the same size.\n";
 		std::cerr<<"                                                : Tree from "<<fileName1_<<" has "<<size1<<" experiments.";
 		std::cerr<<"                                                : Tree from "<<fileName2_<<" has "<<size2<<" experiments.";
 		return kFALSE;
 	}
 
 	for ( ExptsMap::const_iterator iter1 = tree1Expts_.begin(); iter1 != tree1Expts_.end(); ++iter1 ) {
 		Int_t expt = iter1->first;
 		ExptsMap::const_iterator iter2 = tree2Expts_.find( expt );
 		if ( iter2 == tree2Expts_.end() ) {
 			std::cerr<<"ERROR in LauMergeDataFiles::checkExperimentMaps : Cannot find experiment "<<expt<<" in tree from "<<fileName2_<<std::endl;
 			return kFALSE;
 		}
 	}
 
 	return kTRUE;
 }
 
 void LauMergeDataFiles::readExperiment(TTree* tree, const ExptsMap::const_iterator& expt, Int_t offset)
 {
 	// find the first and last entry for this experiment from the map element
 	const Int_t firstEntry = expt->second.first;
 	const Int_t lastEntry = expt->second.second;
 
 	// loop through all the entries
 	for ( Int_t iEntry(firstEntry); iEntry<=lastEntry; ++iEntry ) {
 		// get the entry from the tree
 		tree->GetEntry( iEntry );
 		// apply the offset to the "event within experiment" variable
 		iEvtWithinExpt_ += offset;
 		// fill the output tree
 		outputTree_->Fill();
 	}
 }
 
 void LauMergeDataFiles::writeFile()
 {
 	std::cout<<"Building experiment:event index"<<std::endl;
 	outputTree_->BuildIndex("iExpt","iEvtWithinExpt");
 
 	std::cout<<"Writing data to outputfile "<<outputFile_->GetName()<<std::endl;
 	outputTree_->SetDirectory(outputFile_);
 	outputFile_->Write();
 
 	// clean-up
 	outputFile_->Close();
 	delete outputFile_; outputFile_ = 0; outputTree_ = 0;
 
 	inputFile1_->Close();
 	delete inputFile1_; inputFile1_ = 0; inputTree1_ = 0;
 
 	inputFile2_->Close();
 	delete inputFile2_; inputFile2_ = 0; inputTree2_ = 0;
 
 	doubleVars_.clear();
 	integerVars_.clear();
 	tree1Expts_.clear();
 	tree2Expts_.clear();
 }
 
diff --git a/src/LauResultsExtractor.cc b/src/LauResultsExtractor.cc
index fb46651..f4e424d 100644
--- a/src/LauResultsExtractor.cc
+++ b/src/LauResultsExtractor.cc
@@ -1,292 +1,296 @@
 
 /*
 Copyright 2005 University of Warwick
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at
 
     http://www.apache.org/licenses/LICENSE-2.0
 
 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.
 */
 
 /*
 Laura++ package authors:
 John Back
 Paul Harrison
 Thomas Latham
 */
 
+/*! \file LauResultsExtractor.cc
+  \brief File containing implementation of LauResultsExtractor class.
+ */
+
 #include "LauResultsExtractor.hh"
 
 #include <iostream>
 #include <fstream>
 #include <cstdlib>
 #include <vector>
 #include <map>
 
 #include "TChain.h"
 #include "TFile.h"
 #include "TH1.h"
 #include "TLeaf.h"
 #include "TObjArray.h"
 #include "TSystem.h"
 
 ClassImp(LauResultsExtractor)
 
 
 LauResultsExtractor::LauResultsExtractor(const TString& inputFileName, const TString& outputFileName, const TString& treeName) :
 	inputFileName_(inputFileName),
 	outputFileName_(outputFileName),
 	treeName_(treeName),
 	inputTree_(0),
 	outputFile_(0),
 	outputTree_(0),
 	nEntries_(0)
 {
 }
 
 LauResultsExtractor::~LauResultsExtractor()
 {
 	this->clearMaps();
 	delete inputTree_; inputTree_ = 0;
 	if (outputFile_ && outputFile_->IsOpen()) {
 		delete outputTree_; outputTree_ = 0;
 	}
 	delete outputFile_; outputFile_ = 0;
 }
 
 void LauResultsExtractor::setupInputTree()
 {
 	TObjArray* leaves = inputTree_->GetListOfLeaves();
 	Int_t nLeaves = leaves->GetEntries();
 
 	std::cout<<"Setting branches for input tree \""<<inputTree_->GetName()<<"\" with "<<nEntries_<<" entries..."<<std::endl;
 	inputTree_->SetBranchAddress("iExpt",&iExpt_);
 	inputTree_->SetBranchAddress("fitStatus",&fitStatus_);
 	inputTree_->SetBranchAddress("NLL",&NLL_);
 	inputTree_->SetBranchAddress("EDM",&EDM_);
 
 	for (Int_t iLeaf(3); iLeaf<nLeaves; ++iLeaf) {
 
 		TLeaf * leaf = dynamic_cast<TLeaf*>((*leaves)[iLeaf]);
 		TString type = leaf->GetTypeName();
 		TString name = leaf->GetName();
 		Int_t size = leaf->GetNdata();
 
 		if ((type != "Double_t") || (size != 1)) {
 			continue;
 		}
 
 		std::pair<std::map<TString,Double_t>::iterator,bool> result = otherVars_.insert(std::make_pair(name,0.0));
 		std::map<TString,Double_t>::iterator iter = result.first;
 		bool ok = result.second;
 		if (ok) {
 			inputTree_->SetBranchAddress(name,&(iter->second));
 		}
 	}
 
 	std::cout<<"Set branch addresses for "<<otherVars_.size()+3<<" branches.\n"<<std::endl;
 }
 
 void LauResultsExtractor::setupOutputTree(TTree * tree)
 {
 	std::cout<<"Creating branches for output tree \""<<tree->GetName()<<"\"..."<<std::endl;
 
 	tree->Branch("iExpt",&iExpt_,"iExpt/I");
 	tree->Branch("fitStatus",&fitStatus_,"fitStatus/I");
 	tree->Branch("NLL",&NLL_,"NLL/D");
 	tree->Branch("EDM",&EDM_,"EDM/D");
 
 	for (std::map<TString,Double_t>::iterator iter = otherVars_.begin(); iter != otherVars_.end(); ++iter) {
 		TString name = iter->first;
 		Double_t * address = &(iter->second);
 		TString thirdBit = name;
 		thirdBit += "/D";
 
 		tree->Branch(name,address,thirdBit);
 	}
 	std::cout<<"Created "<<otherVars_.size()<<" branches.\n"<<std::endl;
 }
 
 void LauResultsExtractor::setInputTreeBranchStatus(const Bool_t status)
 {
 	inputTree_->SetBranchStatus("iExpt",kTRUE);
 	inputTree_->SetBranchStatus("fitStatus",kTRUE);
 	inputTree_->SetBranchStatus("NLL",kTRUE);
 	inputTree_->SetBranchStatus("EDM",kTRUE);
 
 	for (std::map<TString,Double_t>::iterator iter = otherVars_.begin(); iter != otherVars_.end(); ++iter) {
 		TString name = iter->first;
 		inputTree_->SetBranchStatus(name,status);
 	}
 }
 
 void LauResultsExtractor::clearMaps()
 {
 	for (std::map<Int_t,TH1*>::iterator iter = nllHistos_.begin(); iter != nllHistos_.end(); ++iter) {
 		delete (iter->second);
 	}
 	bestNLL_.clear();
 	worstNLL_.clear();
 	allNLLs_.clear();
 	nllHistos_.clear();
 }
 
 void LauResultsExtractor::process(Int_t numExpts)
 {
 	// open the text file
 	std::cout << "\n" << "Chaining...\n" << std::endl;
 	std::ifstream textFile(inputFileName_, std::ios::in);
 	if (!textFile.good()) {
 		std::cerr<<"Problem opening file: \""<<inputFileName_<<"\", exiting..."<<std::endl;
 		gSystem->Exit(EXIT_FAILURE);
 	}
 
 	if (inputTree_) { delete inputTree_;  inputTree_ = 0; }
 	inputTree_ = new TChain(treeName_);
 
 	// Read the text file and add each valid entry to the chain
 	TString inputFileName = "";
 	while(inputFileName.ReadLine(textFile) && (!inputFileName.IsNull())) {
 		if (inputFileName.EndsWith(".root") && !inputFileName.BeginsWith("#")) {
 			std::cout << inputFileName << std::endl;
 			inputTree_->Add(inputFileName);
 		}
 		else {
 			std::cout << inputFileName << "\t *** Skipped ***" << std::endl;
 		}
 	}
 
 	textFile.close();
 	std::cout << "\n" << "... finished.\n" << std::endl;
 
 	nEntries_ = inputTree_->GetEntries();
 	this->setupInputTree();
 
 	outputTree_ = new TTree(treeName_,"");
 	this->setupOutputTree(outputTree_);
 
 	// setup the map:
 	// for each experiment there is a pair object holding
 	// the best NLL and the tree entry for that NLL value
 	// each expt starts out with NLL = 0.0 and entry = -1
 	std::cout<<"Setting up the map..."<<std::flush;
 	this->clearMaps();
 	for (Int_t i(0); i<numExpts; ++i) {
 		bestNLL_.insert(std::make_pair(i, std::make_pair(0.0,-1)));
 		worstNLL_.insert(std::make_pair(i, std::make_pair(0.0,-1)));
 		allNLLs_.insert(std::make_pair(i, std::vector<Double_t>()));
 		allNLLs_[i].reserve(nEntries_);
 	}
 	std::cout<<" done.\n"<<std::endl;
 
 	// only read the 3 info branches
 	//this->setInputTreeBranchStatus(kFALSE);
 
 	// loop over the tree and store the best entries for each expt
 	std::cout<<"Starting to store best entry info..."<<std::endl;
 	for (Int_t j(0); j<nEntries_; ++j) {
 
 		if ((nEntries_<100) || (j%(nEntries_/100)==0)) {
 			std::cout<<"Examining entry "<<j<<std::endl;
 		}
 
 		inputTree_->GetEntry(j);
 
 		if ( (fitStatus_ == 3) && (NLL_ > -DBL_MAX/10.0) ) {
 			allNLLs_[iExpt_].push_back(NLL_);
 
 			Double_t curBestNLL = bestNLL_[iExpt_].first;
 			Int_t curBestEntry = bestNLL_[iExpt_].second;
 			if ((NLL_ < curBestNLL) || (curBestEntry == -1)) {
 				bestNLL_[iExpt_] = std::make_pair(NLL_,j);
 			}
 
 			Double_t curWorstNLL = worstNLL_[iExpt_].first;
 			Int_t curWorstEntry = worstNLL_[iExpt_].second;
 			if ((NLL_ > curWorstNLL) || (curWorstEntry == -1)) {
 				worstNLL_[iExpt_] = std::make_pair(NLL_,j);
 			}
 		}
 
 	}
 	std::cout<<"Finished storing best entry info.\n"<<std::endl;
 
 	std::cout<<"Creating NLL histograms..."<<std::flush;
 	TH1* histo(0);
 	for (Int_t i(0); i<numExpts; ++i) {
 		Double_t min = bestNLL_[i].first;
 		Double_t max = worstNLL_[i].first;
 		Double_t range = max - min;
 		if (range < 1e-3) {
 			min -= 0.005;
 			max += 0.005;
 		} else {
 			min -= range*0.2;
 			max += range*0.2;
 		}
 		TString name("expt");
 		name += i;
 		name += "NLL";
 		histo = new TH1F(name,"",100,min,max);
 		for (std::vector<Double_t>::const_iterator iter = allNLLs_[i].begin(); iter != allNLLs_[i].end(); ++iter) {
 			histo->Fill(*iter);
 		}
 		nllHistos_.insert(std::make_pair(i, histo));
 	}
 	std::cout<<" done.\n"<<std::endl;
 
 	// now need to read all branches
 	//this->setInputTreeBranchStatus(kTRUE);
 
 	std::ofstream fout("best-fit.txt");
 
 	// loop over the experiments, grab the best entry and store it
 	std::cout<<"Starting to retrieve best entries and fill output tree."<<std::endl;
 	for (Int_t i(0); i<numExpts; ++i) {
 		Int_t bestEntry = bestNLL_[i].second;
 		if (bestEntry != -1) {
 			inputTree_->GetEntry(bestEntry);
 			outputTree_->Fill();
 		}
 		if ((numExpts<100) || (i%(numExpts/100)==0)) {
 			std::cout<<"Writing out experiment "<<i<<std::endl;
 		}
 		TString bestFit(inputTree_->GetCurrentFile()->GetName());
 		bestFit.Remove(0,3);
 		Int_t index = bestFit.Index("_");
 		if ( index < 1 ) {
 			index = bestFit.Index(".");
 		}
 		bestFit.Remove(index);
 		fout<<"Experiment "<<i<<" BestFit "<<bestFit<<std::endl;
 	}
 	std::cout<<"Finished filling best entries in output tree.\n"<<std::endl;
 	fout.close();
 
 	std::cout<<"Writing output file."<<std::endl;
 
 	this->writeFile();
 }
 
 void LauResultsExtractor::writeFile()
 {
 	if (!outputFile_) {
 		outputFile_ = new TFile(outputFileName_,"recreate");
 	}
 	for (std::map<Int_t,TH1*>::iterator iter = nllHistos_.begin(); iter != nllHistos_.end(); ++iter) {
 		(iter->second)->SetDirectory(outputFile_);
 	}
 	outputTree_->SetDirectory(outputFile_);
 	outputFile_->Write();
 	outputFile_->Close();
 	delete outputFile_; outputFile_ = 0;
 	nllHistos_.clear();
 }