Index: trunk/src/LauResultsExtractor.cc =================================================================== --- trunk/src/LauResultsExtractor.cc (revision 0) +++ trunk/src/LauResultsExtractor.cc (revision 562) @@ -0,0 +1,295 @@ + +/* +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 +*/ + +#include "LauResultsExtractor.hh" + +#include +#include +#include +#include +#include + +#include "TChain.h" +#include "TFile.h" +#include "TH1.h" +#include "TLeaf.h" +#include "TObjArray.h" +#include "TSystem.h" + +using std::cout; +using std::cerr; +using std::flush; +using std::endl; +using std::vector; +using std::map; +using std::pair; +using std::ofstream; + +ResultsExtractor::ResultsExtractor(const TString& inputFileName, const TString& outputFileName, const TString& treeName) : + inputFileName_(inputFileName), + outputFileName_(outputFileName), + treeName_(treeName), + inputTree_(0), + outputFile_(0), + outputTree_(0), + nEntries_(0) +{ +} + +ResultsExtractor::~ResultsExtractor() +{ + this->clearMaps(); + delete inputTree_; inputTree_ = 0; + if (outputFile_ && outputFile_->IsOpen()) { + delete outputTree_; outputTree_ = 0; + } + delete outputFile_; outputFile_ = 0; +} + +void ResultsExtractor::setupInputTree() +{ + TObjArray* leaves = inputTree_->GetListOfLeaves(); + Int_t nLeaves = leaves->GetEntries(); + + cout<<"Setting branches for input tree \""<GetName()<<"\" with "<SetBranchAddress("iExpt",&iExpt_); + inputTree_->SetBranchAddress("fitStatus",&fitStatus_); + inputTree_->SetBranchAddress("NLL",&NLL_); + + for (Int_t iLeaf(3); iLeaf((*leaves)[iLeaf]); + TString type = leaf->GetTypeName(); + TString name = leaf->GetName(); + Int_t size = leaf->GetNdata(); + + if ((type != "Double_t") || (size != 1)) { + continue; + } + + std::pair::iterator,bool> result = otherVars_.insert(std::make_pair(name,0.0)); + std::map::iterator iter = result.first; + bool ok = result.second; + if (ok) { + inputTree_->SetBranchAddress(name,&(iter->second)); + } + } + + cout<<"Set branch addresses for "<GetName()<<"\"..."<Branch("iExpt",&iExpt_,"iExpt/I"); + tree->Branch("fitStatus",&fitStatus_,"fitStatus/I"); + tree->Branch("NLL",&NLL_,"NLL/D"); + + for (std::map::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); + } + cout<<"Created "<SetBranchStatus("iExpt",kTRUE); + inputTree_->SetBranchStatus("fitStatus",kTRUE); + inputTree_->SetBranchStatus("NLL",kTRUE); + + for (std::map::iterator iter = otherVars_.begin(); iter != otherVars_.end(); ++iter) { + TString name = iter->first; + inputTree_->SetBranchStatus(name,status); + } +} + +void ResultsExtractor::clearMaps() +{ + for (map::iterator iter = nllHistos_.begin(); iter != nllHistos_.end(); ++iter) { + delete (iter->second); + } + bestNLL_.clear(); + worstNLL_.clear(); + allNLLs_.clear(); + nllHistos_.clear(); +} + +void ResultsExtractor::process(Int_t numExpts) +{ + // open the text file + cout << "\n" << "Chaining...\n" << endl; + std::ifstream textFile(inputFileName_, std::ios::in); + if (!textFile.good()) { + cerr<<"Problem opening file: \""<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("#")) { + cout << inputFileName << endl; + inputTree_->Add(inputFileName); + } + else { + cout << inputFileName << "\t *** Skipped ***" << endl; + } + } + + textFile.close(); + cout << "\n" << "... finished.\n" << 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 + cout<<"Setting up the map..."<clearMaps(); + for (Int_t i(0); i())); + allNLLs_[i].reserve(nEntries_); + } + cout<<" done.\n"<setInputTreeBranchStatus(kFALSE); + + // loop over the tree and store the best entries for each expt + cout<<"Starting to store best entry info..."<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); + } + } + + } + cout<<"Finished storing best entry info.\n"<::const_iterator iter = allNLLs_[i].begin(); iter != allNLLs_[i].end(); ++iter) { + histo->Fill(*iter); + } + nllHistos_.insert(std::make_pair(i, histo)); + } + cout<<" done.\n"<setInputTreeBranchStatus(kTRUE); + + std::ofstream fout("best-fit.txt"); + + // loop over the experiments, grab the best entry and store it + cout<<"Starting to retrieve best entries and fill output tree."<GetEntry(bestEntry); + outputTree_->Fill(); + } + if ((numExpts<100) || (i%(numExpts/100)==0)) { + cout<<"Writing out experiment "<GetCurrentFile()->GetName()); + bestFit.Remove(0,3); + Int_t index = bestFit.Index("_"); + if ( index < 1 ) { + index = bestFit.Index("."); + } + bestFit.Remove(index); + fout<<"Experiment "<writeFile(); +} + +void ResultsExtractor::writeFile() +{ + if (!outputFile_) { + outputFile_ = new TFile(outputFileName_,"recreate"); + } + for (map::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(); +} + Index: trunk/src/LauCalcChiSq.cc =================================================================== --- trunk/src/LauCalcChiSq.cc (revision 0) +++ trunk/src/LauCalcChiSq.cc (revision 562) @@ -0,0 +1,546 @@ + +/* +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 + +#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 +#include +#include +#include + +using std::cout; +using std::cerr; +using std::endl; + +CalcChiSq::CalcChiSq(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) + +{ +} + +CalcChiSq::~CalcChiSq() +{ +} + +void CalcChiSq::run() +{ + + cout<<"Running chi-squared algorithm"<initialiseHistos(); + cout<<"Calculating chi-squared"<calculateChiSq(); + + //make plots + makePlots(); + + // Output the various histograms + cout<<"Writing out histogram output"<SetDirectory(outFile); + histo2_->SetDirectory(outFile); + pullHisto_->SetDirectory(outFile); + chiSqHisto_->SetDirectory(outFile); + chiSqSignedHisto_->SetDirectory(outFile); + + outFile->Write(); + outFile->Close(); + + cout<<"Done"<GetBinContent(1); + + for(Int_t i=1; i<=histo1_->GetNumberOfBins(); ++i) { + //keep track of actual minimum + if(histo1_->GetBinContent(i)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; + + cout<<"Total ChiSq/nDof = "<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 CalcChiSq::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()) { + 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) { + cerr<<"Error. Could not find the tree "<Exit(EXIT_FAILURE); + } + + // get the info on the high stat histogram + getData >> fileName2_ >> treeName2_ >> xName2_ >> yName2_; + if (!getData.good()) { + cerr<<"Error. Could not read the second line of the input file "<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(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(file2->Get(treeName2_.Data())); + } + + if (tree2 == 0) { + cerr<<"Error. Could not find the tree "<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(); + + cout<<"Using the files and trees: "<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 divisions; + pickBinning(xs,ys,nEntries,divisions); + + //perform the adaptive bining based on histo1_ + getHisto(xMin_, xMax_, yMin_, yMax_, xs, ys, nEntries, 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); + + 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 CalcChiSq::pickBinning(Double_t* xs, Double_t* ys, Int_t nEntries, std::vector& divisions) { + //first check how many events we have within the histogram limits + Int_t nIn(0); + for(Int_t i=0; i= xMin_ && 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& divisions, 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; ixMax)||(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; ixMax)||(ys[i]yMax)) continue; + if ((xs[i]=xLimits[nDivx+1])||(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){ + getHisto(xLimits[nDivx], xLimits[nDivx + 1], yLimits[nDivx][nDivy], yLimits[nDivx][nDivy + 1], xs, ys, nEntries, divisions,iter+1); + } + } +} + Index: trunk/src/LauMergeDataFiles.cc =================================================================== --- trunk/src/LauMergeDataFiles.cc (revision 0) +++ trunk/src/LauMergeDataFiles.cc (revision 562) @@ -0,0 +1,303 @@ + +/* +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 +*/ + +#include "LauMergeDataFiles.hh" + +#include +#include + +#include "TLeaf.h" +#include "TObjArray.h" +#include "TSystem.h" + +MergeDataFiles::MergeDataFiles(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) +{ +} + +MergeDataFiles::~MergeDataFiles() +{ + 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 MergeDataFiles::openInputFiles() +{ + // open the two input ROOT files + inputFile1_ = TFile::Open(fileName1_); + if (!inputFile1_ || inputFile1_->IsZombie()) { + std::cerr<<"Problem opening file: \""<Exit(EXIT_FAILURE); + } + inputTree1_ = dynamic_cast( inputFile1_->Get(treeName_) ); + if (!inputTree1_) { + std::cerr<<"Problem getting tree called "<Exit(EXIT_FAILURE); + } + + inputFile2_ = TFile::Open(fileName2_); + if (!inputFile2_ || inputFile2_->IsZombie()) { + std::cerr<<"Problem opening file: \""<Exit(EXIT_FAILURE); + } + inputTree2_ = dynamic_cast( inputFile2_->Get(treeName_) ); + if (!inputTree2_) { + std::cerr<<"Problem getting tree called "<Exit(EXIT_FAILURE); + } +} + +void MergeDataFiles::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."<SetBranchAddress("iExpt",&iExpt_); + inputTree1_->SetBranchAddress("iEvtWithinExpt",&iEvtWithinExpt_); + inputTree2_->SetBranchAddress("iExpt",&iExpt_); + inputTree2_->SetBranchAddress("iEvtWithinExpt",&iEvtWithinExpt_); + + for (Int_t iLeaf(0); iLeaf((*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 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 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 "<GetName()<<"\"..."<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 "<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..."<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 MergeDataFiles::findExperiments(TTree* tree, ExptsMap& exptsMap) +{ + const Int_t nEntries = tree->GetEntries(); + + // loop through the tree + for ( Int_t iEntry(0); iEntryGetEntry(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 MergeDataFiles::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 MergeDataFiles::checkExperimentMaps : Experiment maps are not the same size.\n"; + std::cerr<<" : Tree from "<first; + ExptsMap::const_iterator iter2 = tree2Expts_.find( expt ); + if ( iter2 == tree2Expts_.end() ) { + std::cerr<<"ERROR in MergeDataFiles::checkExperimentMaps : Cannot find experiment "<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 MergeDataFiles::writeFile() +{ + std::cout<<"Building experiment:event index"<BuildIndex("iExpt","iEvtWithinExpt"); + + std::cout<<"Writing data to outputfile "<GetName()<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(); +} + Index: trunk/inc/LauResultsExtractor.hh =================================================================== --- trunk/inc/LauResultsExtractor.hh (revision 0) +++ trunk/inc/LauResultsExtractor.hh (revision 562) @@ -0,0 +1,85 @@ + +/* +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 +*/ + +#include +#include + +#include "TString.h" + +class TChain; +class TFile; +class TH1; +class TTree; + +/* + A utility class to allow the extraction of the best fit from a series of + fits to a given data sample. + When fitting amplitude models, the likelihood parameter space is highly + complex. Hence the fitter can often wander into local minima. To mitigate + this effect a data sample can be fitted many times with randomised starting + values of the isobar parameters. It is then necessary to determine which of + these fits gives the best solution, i.e. the minimum NLL. This class + performs this task, reading in a series of fits for each data sample and + writing out a single file that contains the results of the best fit for each + data sample. +*/ + +class ResultsExtractor +{ + public: + ResultsExtractor(const TString& inputFileName, const TString& outputFileName, const TString& treeName); + ~ResultsExtractor(); + + void process(Int_t numExpts); + + protected: + void setupInputTree(); + void setupOutputTree(TTree * tree); + void setInputTreeBranchStatus(Bool_t status); + void clearMaps(); + void writeFile(); + + private: + TString inputFileName_; + TString outputFileName_; + TString treeName_; + TChain * inputTree_; + TFile * outputFile_; + TTree * outputTree_; + Int_t nEntries_; + + // Tree variables + Int_t iExpt_; + Int_t fitStatus_; + Double_t NLL_; + + std::map otherVars_; + + std::map< Int_t, std::pair > bestNLL_; + std::map< Int_t, std::pair > worstNLL_; + std::map< Int_t, std::vector > allNLLs_; + + std::map< Int_t, TH1* > nllHistos_; +}; + Index: trunk/inc/LauMergeDataFiles.hh =================================================================== --- trunk/inc/LauMergeDataFiles.hh (revision 0) +++ trunk/inc/LauMergeDataFiles.hh (revision 562) @@ -0,0 +1,81 @@ + +/* +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 +*/ + +#include +#include + +#include "TFile.h" +#include "TString.h" +#include "TTree.h" + +// A utility class to allow the merging of data files on a expt-by-expt basis, +// i.e. such that events for expt 0 from tree 1 will be followed by events for expt 0 from tree 2, +// then expt 1 from tree1, expt 1 from tree 2 etc. + +class MergeDataFiles +{ + public: + MergeDataFiles(const TString& fileName1, const TString& fileName2, const TString& treeName); + ~MergeDataFiles(); + + void process(const TString& fileName); + + protected: + typedef std::map LeafDoubleMap; + typedef std::map LeafIntegerMap; + typedef std::map< Int_t,std::pair > ExptsMap; + + void openInputFiles(); + void setupInputTrees(); + void setupOutputTree(); + void findExperiments(TTree* tree, ExptsMap& exptsMap); + Bool_t checkExperimentMaps() const; + void readExperiment(TTree* tree, const ExptsMap::const_iterator& exptsMap, Int_t offset); + void writeFile(); + + private: + TString fileName1_; + TString fileName2_; + TString treeName_; + + TFile * inputFile1_; + TFile * inputFile2_; + + TTree * inputTree1_; + TTree * inputTree2_; + + TFile * outputFile_; + TTree * outputTree_; + + // Tree variables + Int_t iExpt_; + Int_t iEvtWithinExpt_; + + LeafDoubleMap doubleVars_; + LeafIntegerMap integerVars_; + + ExptsMap tree1Expts_; + ExptsMap tree2Expts_; +}; + Index: trunk/inc/LauCalcChiSq.hh =================================================================== --- trunk/inc/LauCalcChiSq.hh (revision 0) +++ trunk/inc/LauCalcChiSq.hh (revision 562) @@ -0,0 +1,73 @@ + +/* +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 +*/ + +#include "TH2Poly.h" +#include "TString.h" + +#include + +/* + A utility class to allow the calculation of the chisq of the fit to the Dalitz plot. + A text config file is provided that gives the datasets for the data and + toy MC generated from the fit results. These can be in the traditional DP + or the square DP. A sample config file is provided. The fields are: + + + +*/ + +class CalcChiSq { + +public: + + CalcChiSq(TString inputFileName = "chiSqInput.txt"); + virtual ~CalcChiSq(); + + void run(); + inline void setVerbose(Bool_t flag) {verbose_ = flag;} + +private: + + void initialiseHistos(); + void pickBinning(Double_t* xs, Double_t* ys, Int_t nEntries, std::vector& divisions); + void getHisto(Double_t xMin, Double_t xMax, Double_t yMin, Double_t yMax, Double_t* xs, Double_t* ys, Int_t nEntries, std::vector& divisions, UInt_t iter=0); + + void calculateChiSq(); + void makePlots(); + + TString inputFileName_; + TString fileName1_, fileName2_, treeName1_, treeName2_, xName1_, xName2_, yName1_, yName2_; + Float_t minContent_; + + TH2Poly *theHisto_, *histo1_, *histo2_; + TH2Poly *pullHisto_; + TH2Poly *chiSqHisto_; + TH2Poly *chiSqSignedHisto_; + + Float_t xMin_, xMax_, yMin_, yMax_; + Int_t nParams_; + Float_t scaleFactor_; + Bool_t verbose_; + +}; Index: trunk/examples/ResultsExtractorMain.cc =================================================================== --- trunk/examples/ResultsExtractorMain.cc (revision 561) +++ trunk/examples/ResultsExtractorMain.cc (revision 562) @@ -1,55 +0,0 @@ - -/* -Copyright 2013 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 -*/ - -#include - -#include "TString.h" - -#include "ResultsExtractor.hh" - -int main(const int argc, const char** argv) -{ - Int_t nExpt(1); - TString inputFileName("input-list.txt"); - TString outputFileName("output.root"); - TString treeName("fitResults"); - - if (argc > 1) { - nExpt = atoi(argv[1]); - } - if (argc > 2) { - inputFileName = argv[2]; - } - if (argc > 3) { - outputFileName = argv[3]; - } - if (argc > 4) { - treeName = argv[4]; - } - - ResultsExtractor myExtractor(inputFileName,outputFileName,treeName); - myExtractor.process(nExpt); - - return EXIT_SUCCESS; -} Index: trunk/examples/ResultsExtractor.hh =================================================================== --- trunk/examples/ResultsExtractor.hh (revision 561) +++ trunk/examples/ResultsExtractor.hh (revision 562) @@ -1,85 +0,0 @@ - -/* -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 -*/ - -#include -#include - -#include "TString.h" - -class TChain; -class TFile; -class TH1; -class TTree; - -/* - A utility class to allow the extraction of the best fit from a series of - fits to a given data sample. - When fitting amplitude models, the likelihood parameter space is highly - complex. Hence the fitter can often wander into local minima. To mitigate - this effect a data sample can be fitted many times with randomised starting - values of the isobar parameters. It is then necessary to determine which of - these fits gives the best solution, i.e. the minimum NLL. This class - performs this task, reading in a series of fits for each data sample and - writing out a single file that contains the results of the best fit for each - data sample. -*/ - -class ResultsExtractor -{ - public: - ResultsExtractor(const TString& inputFileName, const TString& outputFileName, const TString& treeName); - ~ResultsExtractor(); - - void process(Int_t numExpts); - - protected: - void setupInputTree(); - void setupOutputTree(TTree * tree); - void setInputTreeBranchStatus(Bool_t status); - void clearMaps(); - void writeFile(); - - private: - TString inputFileName_; - TString outputFileName_; - TString treeName_; - TChain * inputTree_; - TFile * outputFile_; - TTree * outputTree_; - Int_t nEntries_; - - // Tree variables - Int_t iExpt_; - Int_t fitStatus_; - Double_t NLL_; - - std::map otherVars_; - - std::map< Int_t, std::pair > bestNLL_; - std::map< Int_t, std::pair > worstNLL_; - std::map< Int_t, std::vector > allNLLs_; - - std::map< Int_t, TH1* > nllHistos_; -}; - Index: trunk/examples/MergeDataFilesMain.cc =================================================================== --- trunk/examples/MergeDataFilesMain.cc (revision 561) +++ trunk/examples/MergeDataFilesMain.cc (revision 562) @@ -1,55 +0,0 @@ - -/* -Copyright 2013 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 -*/ - -#include - -#include "TString.h" - -#include "MergeDataFiles.hh" - -int main(const int argc, const char** argv) -{ - if (argc < 3) { - std::cerr<<"Usage: "< [treeName] [outputFileName]"< 3) { - treeName = argv[3]; - } - TString outputFileName("test.root"); - if (argc > 4) { - outputFileName = argv[4]; - } - - MergeDataFiles myMerger(inputFileName1,inputFileName2,treeName); - myMerger.process(outputFileName); - - return EXIT_SUCCESS; -} - Index: trunk/examples/MergeDataFiles.hh =================================================================== --- trunk/examples/MergeDataFiles.hh (revision 561) +++ trunk/examples/MergeDataFiles.hh (revision 562) @@ -1,81 +0,0 @@ - -/* -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 -*/ - -#include -#include - -#include "TFile.h" -#include "TString.h" -#include "TTree.h" - -// A utility class to allow the merging of data files on a expt-by-expt basis, -// i.e. such that events for expt 0 from tree 1 will be followed by events for expt 0 from tree 2, -// then expt 1 from tree1, expt 1 from tree 2 etc. - -class MergeDataFiles -{ - public: - MergeDataFiles(const TString& fileName1, const TString& fileName2, const TString& treeName); - ~MergeDataFiles(); - - void process(const TString& fileName); - - protected: - typedef std::map LeafDoubleMap; - typedef std::map LeafIntegerMap; - typedef std::map< Int_t,std::pair > ExptsMap; - - void openInputFiles(); - void setupInputTrees(); - void setupOutputTree(); - void findExperiments(TTree* tree, ExptsMap& exptsMap); - Bool_t checkExperimentMaps() const; - void readExperiment(TTree* tree, const ExptsMap::const_iterator& exptsMap, Int_t offset); - void writeFile(); - - private: - TString fileName1_; - TString fileName2_; - TString treeName_; - - TFile * inputFile1_; - TFile * inputFile2_; - - TTree * inputTree1_; - TTree * inputTree2_; - - TFile * outputFile_; - TTree * outputTree_; - - // Tree variables - Int_t iExpt_; - Int_t iEvtWithinExpt_; - - LeafDoubleMap doubleVars_; - LeafIntegerMap integerVars_; - - ExptsMap tree1Expts_; - ExptsMap tree2Expts_; -}; - Index: trunk/examples/CalcChiSqMain.cc =================================================================== --- trunk/examples/CalcChiSqMain.cc (revision 561) +++ trunk/examples/CalcChiSqMain.cc (revision 562) @@ -1,44 +0,0 @@ - -/* -Copyright 2013 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 -*/ - -#include "CalcChiSq.hh" - -#include "TString.h" - -int main(const int argc, const char** argv) -{ - - TString inputFile("chiSqInput.txt"); - - if (argc > 1) {inputFile = TString(argv[1]);} - - CalcChiSq a(inputFile); - - //a.setVerbose(kTRUE); - a.run(); - - return 0; - -} - Index: trunk/examples/CalcChiSq.hh =================================================================== --- trunk/examples/CalcChiSq.hh (revision 561) +++ trunk/examples/CalcChiSq.hh (revision 562) @@ -1,73 +0,0 @@ - -/* -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 -*/ - -#include "TH2Poly.h" -#include "TString.h" - -#include - -/* - A utility class to allow the calculation of the chisq of the fit to the Dalitz plot. - A text config file is provided that gives the datasets for the data and - toy MC generated from the fit results. These can be in the traditional DP - or the square DP. A sample config file is provided. The fields are: - - - -*/ - -class CalcChiSq { - -public: - - CalcChiSq(TString inputFileName = "chiSqInput.txt"); - virtual ~CalcChiSq(); - - void run(); - inline void setVerbose(Bool_t flag) {verbose_ = flag;} - -private: - - void initialiseHistos(); - void pickBinning(Double_t* xs, Double_t* ys, Int_t nEntries, std::vector& divisions); - void getHisto(Double_t xMin, Double_t xMax, Double_t yMin, Double_t yMax, Double_t* xs, Double_t* ys, Int_t nEntries, std::vector& divisions, UInt_t iter=0); - - void calculateChiSq(); - void makePlots(); - - TString inputFileName_; - TString fileName1_, fileName2_, treeName1_, treeName2_, xName1_, xName2_, yName1_, yName2_; - Float_t minContent_; - - TH2Poly *theHisto_, *histo1_, *histo2_; - TH2Poly *pullHisto_; - TH2Poly *chiSqHisto_; - TH2Poly *chiSqSignedHisto_; - - Float_t xMin_, xMax_, yMin_, yMax_; - Int_t nParams_; - Float_t scaleFactor_; - Bool_t verbose_; - -}; Index: trunk/examples/MergeDataFiles.cc =================================================================== --- trunk/examples/MergeDataFiles.cc (revision 561) +++ trunk/examples/MergeDataFiles.cc (revision 562) @@ -1,303 +1,55 @@ /* -Copyright 2008 University of Warwick +Copyright 2013 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 */ -#include "MergeDataFiles.hh" - #include -#include - -#include "TLeaf.h" -#include "TObjArray.h" -#include "TSystem.h" - -MergeDataFiles::MergeDataFiles(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) -{ -} - -MergeDataFiles::~MergeDataFiles() -{ - 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 MergeDataFiles::openInputFiles() -{ - // open the two input ROOT files - inputFile1_ = TFile::Open(fileName1_); - if (!inputFile1_ || inputFile1_->IsZombie()) { - std::cerr<<"Problem opening file: \""<Exit(EXIT_FAILURE); - } - inputTree1_ = dynamic_cast( inputFile1_->Get(treeName_) ); - if (!inputTree1_) { - std::cerr<<"Problem getting tree called "<Exit(EXIT_FAILURE); - } - - inputFile2_ = TFile::Open(fileName2_); - if (!inputFile2_ || inputFile2_->IsZombie()) { - std::cerr<<"Problem opening file: \""<Exit(EXIT_FAILURE); - } - inputTree2_ = dynamic_cast( inputFile2_->Get(treeName_) ); - if (!inputTree2_) { - std::cerr<<"Problem getting tree called "<Exit(EXIT_FAILURE); - } -} - -void MergeDataFiles::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."<SetBranchAddress("iExpt",&iExpt_); - inputTree1_->SetBranchAddress("iEvtWithinExpt",&iEvtWithinExpt_); - inputTree2_->SetBranchAddress("iExpt",&iExpt_); - inputTree2_->SetBranchAddress("iEvtWithinExpt",&iEvtWithinExpt_); - - for (Int_t iLeaf(0); iLeaf((*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 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 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 "<GetName()<<"\"..."<Branch("iExpt",&iExpt_,"iExpt/I"); - outputTree_->Branch("iEvtWithinExpt",&iEvtWithinExpt_,"iEvtWithinExpt/I"); +#include "TString.h" - for (LeafDoubleMap::iterator iter = doubleVars_.begin(); iter != doubleVars_.end(); ++iter) { - TString name = iter->first; - Double_t * address = &(iter->second); - TString thirdBit = name; - thirdBit += "/D"; +#include "LauMergeDataFiles.hh" - 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 "<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..."<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 MergeDataFiles::findExperiments(TTree* tree, ExptsMap& exptsMap) +int main(const int argc, const char** argv) { - const Int_t nEntries = tree->GetEntries(); - - // loop through the tree - for ( Int_t iEntry(0); iEntryGetEntry(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; - } - } + if (argc < 3) { + std::cerr<<"Usage: "< [treeName] [outputFileName]"<first; - ExptsMap::const_iterator iter2 = tree2Expts_.find( expt ); - if ( iter2 == tree2Expts_.end() ) { - std::cerr<<"ERROR in MergeDataFiles::checkExperimentMaps : Cannot find experiment "< 3) { + treeName = argv[3]; } - - return kTRUE; -} - -void MergeDataFiles::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(); + TString outputFileName("test.root"); + if (argc > 4) { + outputFileName = argv[4]; } -} -void MergeDataFiles::writeFile() -{ - std::cout<<"Building experiment:event index"<BuildIndex("iExpt","iEvtWithinExpt"); + MergeDataFiles myMerger(inputFileName1,inputFileName2,treeName); + myMerger.process(outputFileName); - std::cout<<"Writing data to outputfile "<GetName()<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(); + return EXIT_SUCCESS; } Index: trunk/examples/Makefile =================================================================== --- trunk/examples/Makefile (revision 561) +++ trunk/examples/Makefile (revision 562) @@ -1,167 +1,150 @@ ############################################################################ # 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 # # # ############################################################################ # # # ---------------------------------------- # # Standalone Makefile for Laura++ examples # # ---------------------------------------- # # # # Instructions # # - Review 'external configuration' section below # # to match systems compilers setup # # # # - Make sure the ROOTSYS environment variable is set and points # # to your ROOT release or the root-config script is in your PATH # # # # - run 'make ' # # # # Build targets # # bin - make all examples (default) # # - make the specific example # # clean - delete all intermediate and final build objects # # # ############################################################################ # --- External configuration ---------------------------------- # first check that ROOTSYS is defined ifndef ROOTSYS ROOTSYS := $(shell root-config --prefix) ROOTBINDIR := $(shell root-config --bindir) ifeq ($(ROOTSYS), ) $(error running of root-config failed or reported null value) endif else ROOTBINDIR := $(ROOTSYS)/bin endif # By default, don't build the LauRooFitSlave class since it depends on RooFit # and we don't want to pull in that library if we don't have to. # If you want to build it, just set the SKIPLIST variable to be empty. SKIPLIST = SlaveRooFit.cc ROOTCONFIG := $(ROOTBINDIR)/root-config ARCH := $(shell $(ROOTCONFIG) --arch) PLATFORM := $(shell $(ROOTCONFIG) --platform) INCLUDES = WORKDIR = tmp ifeq ($(findstring linux, $(ARCH)),linux) # This set here should work for Linux. CXX = g++ LD = g++ CXXFLAGS = -g -O2 -Wall -Wextra -Wshadow -Woverloaded-virtual -Werror -fPIC MFLAGS = -MM LDFLAGS = -g SOFLAGS = -shared endif ifeq ($(ARCH),macosx64) # For Mac OS X you may need to put -m64 in CXXFLAGS and SOFLAGS. CXX = g++ LD = g++ CXXFLAGS = -g -O3 -Wall -Wextra -Wshadow -Woverloaded-virtual -Werror -fPIC -m64 MFLAGS = -MM LDFLAGS = -g SOFLAGS = -m64 -dynamiclib -single_module -undefined dynamic_lookup endif # --- Internal configuration ---------------------------------- INCDIR=../inc DEPDIR=$(WORKDIR)/dependencies OBJDIR=$(WORKDIR)/objects ROOTCFLAGS := $(shell $(ROOTCONFIG) --cflags) ROOTLIBS := $(shell $(ROOTCONFIG) --libs) ROOTLIBS += -lEG ROOTLIBS += -lMinuit ROOTLIBS += -lTreePlayer ifeq ($(strip $(SKIPLIST)),) ROOTLIBS += -lRooFitCore ROOTLIBS += -lRooFit endif LAURALIBDIR=$(shell pwd | xargs dirname)/lib LAURALIB = $(LAURALIBDIR)/libLaura++.so -UTILSLIB = $(PWD)/$(WORKDIR)/libLauraUtils.so INCLUDES += -I$(INCDIR) CXXFLAGS += $(INCLUDES) CXXFLAGS += $(ROOTCFLAGS) default: bin # List of all source files CCLIST:=$(filter-out $(SKIPLIST), $(wildcard *.cc)) -# List of all source files that contain main functions -BINCCLIST:=$(shell egrep -l "^[[:space:]]*int[[:space:]]*main\>" $(CCLIST)) - -# List of all source files to be compiled into the library -LIBCCLIST:=$(filter-out $(BINCCLIST), $(CCLIST)) - # List of all object files to build OLIST:=$(patsubst %.cc,%.o,$(addprefix $(OBJDIR)/,$(notdir $(CCLIST)))) -# List of all object files to be combined into library -LIBOLIST:=$(patsubst %.cc,%.o,$(addprefix $(OBJDIR)/,$(notdir $(LIBCCLIST)))) - # List of all dependency files to make DLIST:=$(patsubst %.cc,%.d,$(addprefix $(DEPDIR)/,$(notdir $(CCLIST)))) # List of all binary files to make -BINLIST:=$(patsubst %.cc,%,$(notdir $(BINCCLIST))) +BINLIST:=$(patsubst %.cc,%,$(notdir $(CCLIST))) # Implicit rule making all dependency Makefiles included at the end of this makefile $(DEPDIR)/%.d: %.cc @echo "Making $@" @mkdir -p $(DEPDIR) @set -e; $(CXX) $(MFLAGS) $(CXXFLAGS) $< \ | sed 's#\($(notdir $*)\)\.o[ :]*#$(OBJDIR)/\1.o $@ : #g' > $@; \ [ -s $@ ] || rm -f $@ # Implicit rule to compile all sources $(OBJDIR)/%.o : %.cc @echo "Compiling $<" @mkdir -p $(OBJDIR) @$(CXX) $(CXXFLAGS) -c $< -o $@ -# Rule to combine objects into a shared library -$(UTILSLIB): $(LIBOLIST) - @echo "Making $(UTILSLIB)" - @rm -f $(UTILSLIB) - @$(CXX) $(LIBOLIST) $(SOFLAGS) -o $(UTILSLIB) - # Rule to compile all binaries -% : $(OBJDIR)/%.o $(LAURALIB) $(UTILSLIB) +% : $(OBJDIR)/%.o $(LAURALIB) @echo "Linking $@" - @$(CXX) $(LDFLAGS) $< -o $@ $(LAURALIB) $(UTILSLIB) $(ROOTLIBS) + @$(CXX) $(LDFLAGS) $< -o $@ $(LAURALIB) $(ROOTLIBS) bin: $(BINLIST) clean: rm -f $(BINLIST) - rm -f $(UTILSLIB) rm -rf $(WORKDIR) .PHONY : bin default clean -include $(DLIST) Index: trunk/examples/CalcChiSq.cc =================================================================== --- trunk/examples/CalcChiSq.cc (revision 561) +++ trunk/examples/CalcChiSq.cc (revision 562) @@ -1,546 +1,44 @@ /* -Copyright 2008 University of Warwick +Copyright 2013 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 - -#include "CalcChiSq.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 -#include -#include -#include - -using std::cout; -using std::cerr; -using std::endl; - -CalcChiSq::CalcChiSq(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) +#include "LauCalcChiSq.hh" -{ -} - -CalcChiSq::~CalcChiSq() -{ -} - -void CalcChiSq::run() -{ - - cout<<"Running chi-squared algorithm"<initialiseHistos(); - cout<<"Calculating chi-squared"<calculateChiSq(); - - //make plots - makePlots(); - - // Output the various histograms - cout<<"Writing out histogram output"<SetDirectory(outFile); - histo2_->SetDirectory(outFile); - pullHisto_->SetDirectory(outFile); - chiSqHisto_->SetDirectory(outFile); - chiSqSignedHisto_->SetDirectory(outFile); - - outFile->Write(); - outFile->Close(); +#include "TString.h" - cout<<"Done"<GetBinContent(1); - - for(Int_t i=1; i<=histo1_->GetNumberOfBins(); ++i) { - //keep track of actual minimum - if(histo1_->GetBinContent(i)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)); - } - } + TString inputFile("chiSqInput.txt"); - ndof = histo1_->GetNumberOfBins()-nParams_-1; + if (argc > 1) {inputFile = TString(argv[1]);} - cout<<"Total ChiSq/nDof = "<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"); - -} + CalcChiSq a(inputFile); -void CalcChiSq::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()) { - 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) { - cerr<<"Error. Could not find the tree "<Exit(EXIT_FAILURE); - } - - // get the info on the high stat histogram - getData >> fileName2_ >> treeName2_ >> xName2_ >> yName2_; - if (!getData.good()) { - cerr<<"Error. Could not read the second line of the input file "<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(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(file2->Get(treeName2_.Data())); - } - - if (tree2 == 0) { - cerr<<"Error. Could not find the tree "<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(); - - cout<<"Using the files and trees: "<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 divisions; - pickBinning(xs,ys,nEntries,divisions); - - //perform the adaptive bining based on histo1_ - getHisto(xMin_, xMax_, yMin_, yMax_, xs, ys, nEntries, 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); - - 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 CalcChiSq::pickBinning(Double_t* xs, Double_t* ys, Int_t nEntries, std::vector& divisions) { - //first check how many events we have within the histogram limits - Int_t nIn(0); - for(Int_t i=0; i= xMin_ && 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& divisions, UInt_t iter) { + return 0; - //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; ixMax)||(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; ixMax)||(ys[i]yMax)) continue; - if ((xs[i]=xLimits[nDivx+1])||(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){ - getHisto(xLimits[nDivx], xLimits[nDivx + 1], yLimits[nDivx][nDivy], yLimits[nDivx][nDivy + 1], xs, ys, nEntries, divisions,iter+1); - } - } } Index: trunk/examples/ResultsExtractor.cc =================================================================== --- trunk/examples/ResultsExtractor.cc (revision 561) +++ trunk/examples/ResultsExtractor.cc (revision 562) @@ -1,295 +1,55 @@ /* -Copyright 2005 University of Warwick +Copyright 2013 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 */ -#include "ResultsExtractor.hh" - -#include -#include #include -#include -#include - -#include "TChain.h" -#include "TFile.h" -#include "TH1.h" -#include "TLeaf.h" -#include "TObjArray.h" -#include "TSystem.h" - -using std::cout; -using std::cerr; -using std::flush; -using std::endl; -using std::vector; -using std::map; -using std::pair; -using std::ofstream; - -ResultsExtractor::ResultsExtractor(const TString& inputFileName, const TString& outputFileName, const TString& treeName) : - inputFileName_(inputFileName), - outputFileName_(outputFileName), - treeName_(treeName), - inputTree_(0), - outputFile_(0), - outputTree_(0), - nEntries_(0) -{ -} - -ResultsExtractor::~ResultsExtractor() -{ - this->clearMaps(); - delete inputTree_; inputTree_ = 0; - if (outputFile_ && outputFile_->IsOpen()) { - delete outputTree_; outputTree_ = 0; - } - delete outputFile_; outputFile_ = 0; -} - -void ResultsExtractor::setupInputTree() -{ - TObjArray* leaves = inputTree_->GetListOfLeaves(); - Int_t nLeaves = leaves->GetEntries(); - - cout<<"Setting branches for input tree \""<GetName()<<"\" with "<SetBranchAddress("iExpt",&iExpt_); - inputTree_->SetBranchAddress("fitStatus",&fitStatus_); - inputTree_->SetBranchAddress("NLL",&NLL_); - - for (Int_t iLeaf(3); iLeaf((*leaves)[iLeaf]); - TString type = leaf->GetTypeName(); - TString name = leaf->GetName(); - Int_t size = leaf->GetNdata(); - - if ((type != "Double_t") || (size != 1)) { - continue; - } - - std::pair::iterator,bool> result = otherVars_.insert(std::make_pair(name,0.0)); - std::map::iterator iter = result.first; - bool ok = result.second; - if (ok) { - inputTree_->SetBranchAddress(name,&(iter->second)); - } - } - cout<<"Set branch addresses for "<GetName()<<"\"..."<Branch("iExpt",&iExpt_,"iExpt/I"); - tree->Branch("fitStatus",&fitStatus_,"fitStatus/I"); - tree->Branch("NLL",&NLL_,"NLL/D"); - - for (std::map::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); - } - cout<<"Created "<SetBranchStatus("iExpt",kTRUE); - inputTree_->SetBranchStatus("fitStatus",kTRUE); - inputTree_->SetBranchStatus("NLL",kTRUE); - - for (std::map::iterator iter = otherVars_.begin(); iter != otherVars_.end(); ++iter) { - TString name = iter->first; - inputTree_->SetBranchStatus(name,status); - } -} +#include "TString.h" -void ResultsExtractor::clearMaps() -{ - for (map::iterator iter = nllHistos_.begin(); iter != nllHistos_.end(); ++iter) { - delete (iter->second); - } - bestNLL_.clear(); - worstNLL_.clear(); - allNLLs_.clear(); - nllHistos_.clear(); -} +#include "LauResultsExtractor.hh" -void ResultsExtractor::process(Int_t numExpts) +int main(const int argc, const char** argv) { - // open the text file - cout << "\n" << "Chaining...\n" << endl; - std::ifstream textFile(inputFileName_, std::ios::in); - if (!textFile.good()) { - cerr<<"Problem opening file: \""<Exit(EXIT_FAILURE); - } - - if (inputTree_) { delete inputTree_; inputTree_ = 0; } - inputTree_ = new TChain(treeName_); + Int_t nExpt(1); + TString inputFileName("input-list.txt"); + TString outputFileName("output.root"); + TString treeName("fitResults"); - // 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("#")) { - cout << inputFileName << endl; - inputTree_->Add(inputFileName); - } - else { - cout << inputFileName << "\t *** Skipped ***" << endl; - } + if (argc > 1) { + nExpt = atoi(argv[1]); } - - textFile.close(); - cout << "\n" << "... finished.\n" << 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 - cout<<"Setting up the map..."<clearMaps(); - for (Int_t i(0); i())); - allNLLs_[i].reserve(nEntries_); + if (argc > 2) { + inputFileName = argv[2]; } - cout<<" done.\n"<setInputTreeBranchStatus(kFALSE); - - // loop over the tree and store the best entries for each expt - cout<<"Starting to store best entry info..."<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); - } - } - + if (argc > 3) { + outputFileName = argv[3]; } - cout<<"Finished storing best entry info.\n"<::const_iterator iter = allNLLs_[i].begin(); iter != allNLLs_[i].end(); ++iter) { - histo->Fill(*iter); - } - nllHistos_.insert(std::make_pair(i, histo)); - } - cout<<" done.\n"<setInputTreeBranchStatus(kTRUE); - - std::ofstream fout("best-fit.txt"); - - // loop over the experiments, grab the best entry and store it - cout<<"Starting to retrieve best entries and fill output tree."<GetEntry(bestEntry); - outputTree_->Fill(); - } - if ((numExpts<100) || (i%(numExpts/100)==0)) { - cout<<"Writing out experiment "<GetCurrentFile()->GetName()); - bestFit.Remove(0,3); - Int_t index = bestFit.Index("_"); - if ( index < 1 ) { - index = bestFit.Index("."); - } - bestFit.Remove(index); - fout<<"Experiment "< 4) { + treeName = argv[4]; } - cout<<"Finished filling best entries in output tree.\n"<writeFile(); + return EXIT_SUCCESS; } - -void ResultsExtractor::writeFile() -{ - if (!outputFile_) { - outputFile_ = new TFile(outputFileName_,"recreate"); - } - for (map::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(); -} -