diff --git a/examples/CalcChiSq.cc b/examples/CalcChiSq.cc index 20e8a5e..9823234 100644 --- a/examples/CalcChiSq.cc +++ b/examples/CalcChiSq.cc @@ -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 +#include "LauCalcChiSq.hh" -// Note that the low and high stat histograms must have the same bin axes -// ranges and number of bins. +#include "TString.h" -// 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) - -{ -} - -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() +int main(const int argc, const char** argv) { - // 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; + TString inputFile("chiSqInput.txt"); - 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; - } + if (argc > 1) {inputFile = TString(argv[1]);} - TMath::Sort(yCountIn, yIn, yIndex, false); + CalcChiSq a(inputFile); - for (Int_t nDivy = 1; nDivy < n_divy; ++nDivy){ - yLimits[nDivx][nDivy] = yIn[yIndex[yCountIn*nDivy/n_divy]]; - } - } + //a.setVerbose(kTRUE); + a.run(); - delete[] xIn; - delete[] yIn; - delete[] xIndex; - delete[] yIndex; + return 0; - //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); - } - } } diff --git a/examples/CalcChiSqMain.cc b/examples/CalcChiSqMain.cc deleted file mode 100644 index 468a0fb..0000000 --- a/examples/CalcChiSqMain.cc +++ /dev/null @@ -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; - -} - diff --git a/examples/Makefile b/examples/Makefile index 1154d81..9c75906 100644 --- a/examples/Makefile +++ b/examples/Makefile @@ -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) diff --git a/examples/MergeDataFiles.cc b/examples/MergeDataFiles.cc index 221c509..29e6306 100644 --- a/examples/MergeDataFiles.cc +++ b/examples/MergeDataFiles.cc @@ -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"); - - 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"; +#include "TString.h" - 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; + if (argc < 3) { + std::cerr<<"Usage: "< [treeName] [outputFileName]"<first; - ExptsMap::const_iterator iter2 = tree2Expts_.find( expt ); + TString inputFileName1(argv[1]); + TString inputFileName2(argv[2]); - // 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 ); + TString treeName("fitTree"); + if (argc > 3) { + treeName = argv[3]; } - - // 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; - } - } + TString outputFileName("test.root"); + if (argc > 4) { + outputFileName = argv[4]; } - // 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; + MergeDataFiles myMerger(inputFileName1,inputFileName2,treeName); + myMerger.process(outputFileName); - doubleVars_.clear(); - integerVars_.clear(); - tree1Expts_.clear(); - tree2Expts_.clear(); + return EXIT_SUCCESS; } diff --git a/examples/MergeDataFilesMain.cc b/examples/MergeDataFilesMain.cc deleted file mode 100644 index c1fc99d..0000000 --- a/examples/MergeDataFilesMain.cc +++ /dev/null @@ -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; -} - diff --git a/examples/ResultsExtractor.cc b/examples/ResultsExtractor.cc index 8ac0621..5be7d4f 100644 --- a/examples/ResultsExtractor.cc +++ b/examples/ResultsExtractor.cc @@ -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; - } +#include "TString.h" - 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"); + Int_t nExpt(1); + TString inputFileName("input-list.txt"); + TString outputFileName("output.root"); + TString treeName("fitResults"); - 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); + if (argc > 1) { + nExpt = atoi(argv[1]); } - 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); + if (argc > 2) { + inputFileName = argv[2]; } -} - -void ResultsExtractor::clearMaps() -{ - for (map::iterator iter = nllHistos_.begin(); iter != nllHistos_.end(); ++iter) { - delete (iter->second); + if (argc > 3) { + outputFileName = argv[3]; } - 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 (argc > 4) { + treeName = argv[4]; } - if (inputTree_) { delete inputTree_; inputTree_ = 0; } - inputTree_ = new TChain(treeName_); + ResultsExtractor myExtractor(inputFileName,outputFileName,treeName); + myExtractor.process(nExpt); - // 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(); + 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(); -} - diff --git a/examples/ResultsExtractorMain.cc b/examples/ResultsExtractorMain.cc deleted file mode 100644 index b6c83f3..0000000 --- a/examples/ResultsExtractorMain.cc +++ /dev/null @@ -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; -} diff --git a/examples/CalcChiSq.hh b/inc/LauCalcChiSq.hh similarity index 100% rename from examples/CalcChiSq.hh rename to inc/LauCalcChiSq.hh diff --git a/examples/MergeDataFiles.hh b/inc/LauMergeDataFiles.hh similarity index 100% rename from examples/MergeDataFiles.hh rename to inc/LauMergeDataFiles.hh diff --git a/examples/ResultsExtractor.hh b/inc/LauResultsExtractor.hh similarity index 100% rename from examples/ResultsExtractor.hh rename to inc/LauResultsExtractor.hh diff --git a/examples/CalcChiSq.cc b/src/LauCalcChiSq.cc similarity index 99% copy from examples/CalcChiSq.cc copy to src/LauCalcChiSq.cc index 20e8a5e..c37f604 100644 --- a/examples/CalcChiSq.cc +++ b/src/LauCalcChiSq.cc @@ -1,546 +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 "CalcChiSq.hh" +#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); } } } diff --git a/examples/MergeDataFiles.cc b/src/LauMergeDataFiles.cc similarity index 99% copy from examples/MergeDataFiles.cc copy to src/LauMergeDataFiles.cc index 221c509..f757d9f 100644 --- a/examples/MergeDataFiles.cc +++ b/src/LauMergeDataFiles.cc @@ -1,303 +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 "MergeDataFiles.hh" +#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(); } diff --git a/examples/ResultsExtractor.cc b/src/LauResultsExtractor.cc similarity index 99% copy from examples/ResultsExtractor.cc copy to src/LauResultsExtractor.cc index 8ac0621..779e98c 100644 --- a/examples/ResultsExtractor.cc +++ b/src/LauResultsExtractor.cc @@ -1,295 +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 "ResultsExtractor.hh" +#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(); }