Log-likelihood scan across two parameters. More...

#include <Surface.h>

Inheritance diagram for ana::Surface:

Public Member Functions
	Surface (const IExperiment expt, osc::IOscCalcAdjustable calc, const IFitVar xvar, int nbinsx, double xmin, double xmax, const IFitVar yvar, int nbinsy, double ymin, double ymax, const std::vector< const IFitVar * > &profVars={}, const std::vector< const ISyst * > &profSysts={}, const std::map< const IFitVar *, std::vector< double >> &seedPts={}, const std::vector< SystShifts > &systSeedPts={}, bool parallel=false, Fitter::Precision prec=Fitter::kNormal)

	Surface (const IExperiment expt, osc::IOscCalcAdjustable calc, const FitAxis &xax, const FitAxis &yax, const std::vector< const IFitVar * > &profVars={}, const std::vector< const ISyst * > &profSysts={}, const std::map< const IFitVar *, std::vector< double >> &seedPts={}, const std::vector< SystShifts > &systSeedPts={}, bool parallel=false, Fitter::Precision prec=Fitter::kNormal)

void	Draw () const
	Draw the surface itself. More...

void	DrawBestFit (Color_t color, Int_t marker=kFullCircle) const
	Draw the best fit point. More...

double	MinChi () const

double	GetMinX () const

double	GetMinY () const

void	DrawContour (TH2 *fc, Style_t style, Color_t color, double minchi=-1)

TH2F *	ToTH2 (double minchi=-1) const

void	SetTitle (const char *str)

std::vector< TH2 * >	GetProfiledHists ()
	Maps of the values taken on by the profiled parameters. More...

std::vector< TH2 * >	GetMarginalizedHists ()
	Deprecated. Retained for backwards compatibility. More...

std::vector< TGraph * >	GetGraphs (TH2 *fc, double minchi=-1)
	For expert use, custom painting of contours. More...

void	SaveTo (TDirectory *dir) const

Static Public Member Functions
static std::unique_ptr< Surface >	LoadFrom (TDirectory *dir)

static std::unique_ptr< Surface >	LoadFromMulti (const std::vector< TFile * > &files, const std::string &label)

static std::unique_ptr< Surface >	LoadFromMulti (const std::string &wildcard, const std::string &label)

Protected Member Functions
	Surface ()

void	EnsureAxes () const

void	CheckMask (const std::string &func) const

void	FillSurface (const std::string &progTitle, const IExperiment expt, osc::IOscCalcAdjustable calc, const FitAxis &xax, const FitAxis &yax, const std::vector< const IFitVar * > &profVars, const std::vector< const ISyst * > &profSysts, const std::map< const IFitVar *, std::vector< double >> &seedPts, const std::vector< SystShifts > &systSeedPts)

void	FillSurfacePoint (const IExperiment expt, osc::IOscCalcAdjustable calc, const FitAxis &xax, double x, const FitAxis &yax, double y, const std::vector< const IFitVar * > &profVars, const std::vector< const ISyst * > &profSysts, const std::map< const IFitVar *, std::vector< double >> &seedPts, const std::vector< SystShifts > &systSeedPts)

Protected Attributes
bool	fParallel

Fitter::Precision	fPrec

double	fMinChi

double	fMinX

double	fMinY

TH2F *	fHist

bool	fLogX

bool	fLogY

std::vector< TH2 * >	fProfHists

std::vector< double >	fSeedValues

std::vector< int >	fBinMask

Friends
class	MedianSurface

Detailed Description

Log-likelihood scan across two parameters.

Definition at line 26 of file sbnana/sbnana/CAFAna/Analysis/Surface.h.

Constructor & Destructor Documentation

ana::Surface::Surface	(	const IExperiment *	expt,
		osc::IOscCalcAdjustable *	calc,
		const IFitVar *	xvar,
		int	nbinsx,
		double	xmin,
		double	xmax,
		const IFitVar *	yvar,
		int	nbinsy,
		double	ymin,
		double	ymax,
		const std::vector< const IFitVar * > &	profVars = `{}`,
		const std::vector< const ISyst * > &	profSysts = `{}`,
		const std::map< const IFitVar *, std::vector< double >> &	seedPts = `{}`,
		const std::vector< SystShifts > &	systSeedPts = `{}`,
		bool	parallel = `false`,
		Fitter::Precision	prec = `Fitter::kNormal`
	)

Parameters

expt	The experiment object to draw $\chi^2$ values from
calc	Values for oscillation parameters to be held fixed
xvar	Oscillation parameter to place on the x axis
nbinsx	Number of bins along x axis
xmin	Minimum value of x axis
xmax	Maximum value of x axis
nbinsy	Number of bins along y axis
ymin	Minimum value of y axis
ymax	Maximum value of y axis
profVars	Oscillation parameters to profile over
profSysts	Systematic parameters to profile over
seedPts	Try all combinations of these params as seeds
systSeedPts	Try all of these systematic combinations as seeds
parallel	Use all the cores on this machine? Be careful...

Definition at line 127 of file sbnana/sbnana/CAFAna/Analysis/Surface.cxx.

     : Surface(expt, calc,
               FitAxis(xvar, nbinsx, xmin, xmax),
               FitAxis(yvar, nbinsy, ymin, ymax),
               profVars, profSysts, seedPts, systSeedPts, parallel, prec)
   {
   }

ana::Surface::Surface	(	const IExperiment *	expt,
		osc::IOscCalcAdjustable *	calc,
		const FitAxis &	xax,
		const FitAxis &	yax,
		const std::vector< const IFitVar * > &	profVars = `{}`,
		const std::vector< const ISyst * > &	profSysts = `{}`,
		const std::map< const IFitVar *, std::vector< double >> &	seedPts = `{}`,
		const std::vector< SystShifts > &	systSeedPts = `{}`,
		bool	parallel = `false`,
		Fitter::Precision	prec = `Fitter::kNormal`
	)

Parameters

expt	The experiment object to draw $\chi^2$ values from
calc	Values for oscillation parameters to be held fixed
xax	A FitAxis for the first variable
yax	A FitAxis for the second variable
profVars	Oscillation parameters to profile over
profSysts	Systematic parameters to profile over
seedPts	Try all combinations of these params as seeds
systSeedPts	Try all of these systematic combinations as seeds
parallel	Use all the cores on this machine? Be careful...

Definition at line 38 of file sbnana/sbnana/CAFAna/Analysis/Surface.cxx.

     : fParallel(parallel), fPrec(prec), fLogX(xax.islog), fLogY(yax.islog)
   {
     fHist = ExpandedHistogram(";"+xax.var->LatexName()+";"+yax.var->LatexName(),
                               xax.nbins, xax.min, xax.max, xax.islog,
                               yax.nbins, yax.min, yax.max, yax.islog);
 
     for(unsigned int i = 0; i < profVars.size()+profSysts.size(); ++i){
       std::string title;
       if(i < profVars.size())
         title = profVars[i]->LatexName();
       else
         title = profSysts[i-profVars.size()]->LatexName();
 
       fProfHists.push_back(ExpandedHistogram
                            (title+";"+xax.var->LatexName()+";"+yax.var->LatexName(),
                             xax.nbins, xax.min, xax.max, xax.islog,
                             yax.nbins, yax.min, yax.max, yax.islog));
     }
 
     std::string progTitle = "Filling surface for "+yax.var->ShortName()+" vs "+xax.var->ShortName();
 
     if(!profVars.empty() || !profSysts.empty()){
       progTitle += " (profiling ";
 
       for(const IFitVar* v: profVars) {
         progTitle += v->ShortName() + ", ";
         fSeedValues.push_back( v->GetValue( calc ) );
       }
       for(const ISyst* s: profSysts)
         progTitle += s->ShortName() + ", ";
 
       // Always have one superfluous ", " at the end
       progTitle.resize(progTitle.size()-2);
       progTitle += ")";
     }
 
     FillSurface(progTitle, expt, calc, xax, yax, profVars, profSysts, seedPts, systSeedPts);
 
     // Location of the best minimum found from filled surface
     double minchi = 1e10;
     int minx = xax.nbins/2;
     int miny = yax.nbins/2;
     for(int x = 1; x < xax.nbins+1; ++x){
       for(int y = 1; y < yax.nbins+1; ++y){
         const int bin = (y-1)*fHist->GetNbinsX() + (x-1);
         if(RunningOnGrid() && NumJobs() > 1
            && !std::count(fBinMask.begin(), fBinMask.end(), bin)) continue;
         const double chi = fHist->GetBinContent(x, y);
         if(chi < minchi && !isnan(chi) && !isinf(chi)){
           minchi = chi;
           minx = x;
           miny = y;
         }
       }
     }
 
     std::vector<const IFitVar*> allVars = {xax.var, yax.var};
     allVars.insert(allVars.end(), profVars.begin(), profVars.end());
     Fitter fit(expt, allVars, profSysts);
     fit.SetPrecision(fPrec);
     // Seed from best grid point
     xax.var->SetValue(calc, BinCenter(fHist->GetXaxis(), minx, xax.islog));
     yax.var->SetValue(calc, BinCenter(fHist->GetYaxis(), miny, yax.islog));
     for(int i = 0; i < (int)fSeedValues.size(); ++i) profVars[i]->SetValue( calc, fSeedValues[i] );
     SystShifts systSeed = SystShifts::Nominal();
     fMinChi = fit.Fit(calc, systSeed, seedPts);
     fMinX = xax.var->GetValue(calc);
     fMinY = yax.var->GetValue(calc);
 
     for(int x = 0; x < xax.nbins+2; ++x){
       for(int y = 0; y < yax.nbins+2; ++y){
         fHist->SetBinContent(x, y, fHist->GetBinContent(x, y)-fMinChi);
       }
     }
 
     fHist->SetMinimum(0);
   }

ana::Surface::Surface ( )

inlineprotected

Definition at line 109 of file sbnana/sbnana/CAFAna/Analysis/Surface.h.

109 {};

Member Function Documentation

void ana::Surface::CheckMask ( const std::string & func ) const

protected

Definition at line 718 of file sbnana/sbnana/CAFAna/Analysis/Surface.cxx.

   {
     if(!fBinMask.empty()) {
       std::cout << "Cannot call " << func << "() on a partial surface!" << std::endl;
       abort();
     }
   }

void ana::Surface::Draw ( ) const

Draw the surface itself.

Definition at line 397 of file sbnana/sbnana/CAFAna/Analysis/Surface.cxx.

   {
     EnsureAxes();
 
     fHist->Draw("colz same");
 
     // colz obliterated them
     gPad->RedrawAxis();
 
     gPad->Update();
   }

void ana::Surface::DrawBestFit	(	Color_t	color,
		Int_t	marker = `kFullCircle`
	)		const

Draw the best fit point.

Definition at line 410 of file sbnana/sbnana/CAFAna/Analysis/Surface.cxx.

   {
     CheckMask("DrawBestFit");
     EnsureAxes();
 
     TMarker* mark = new TMarker(fMinX, fMinY, marker);
     mark->SetMarkerSize(1.5);
     mark->SetMarkerColor(color);
     mark->Draw();
     gPad->Update();
   }

void ana::Surface::DrawContour	(	TH2 *	fc,
		Style_t	style,
		Color_t	color,
		double	minchi = `-1`
	)

Parameters

fc	Surface to compare against for this significance level
style	Line style for TAttLine, solid, dotted, dashed etc
color	Line color for TAttLine
minchi	$\chi^2$ of best fit to compare against. Default: best fit from this surface.

Definition at line 480 of file sbnana/sbnana/CAFAna/Analysis/Surface.cxx.

   {
     CheckMask("DrawContour");
     EnsureAxes();
 
     std::vector<TGraph*> gs = GetGraphs(fc, minchi);
 
     for(TGraph* g: gs){
       g->SetLineWidth(3);//2);
       g->SetLineStyle(style);
       g->SetLineColor(color);
       g->Draw("l");
     }
 
     gPad->Update();
   }

void ana::Surface::EnsureAxes ( ) const

protected

Definition at line 342 of file sbnana/sbnana/CAFAna/Analysis/Surface.cxx.

   {
     // Could have a file temporarily open
     DontAddDirectory guard;
 
     // If this pad has already been drawn in, already has axes
     if(gPad && !gPad->GetListOfPrimitives()->IsEmpty()) return;
 
     // Old, hackier solution
     /*
     std::cout << gPad->GetListOfPrimitives()->GetEntries() << std::endl;
     // Which pads have we already drawn axes in? Never draw axes in them
     // again. Unfortunately UniqueID() never seems to be set. If that's the
     // case, set it to a random value and hope...
     static std::set<UInt_t> already;
     if(already.count(gPad->GetUniqueID())) return;
     if(gPad->GetUniqueID() == 0) gPad->SetUniqueID(rand());
     already.insert(gPad->GetUniqueID());
     */
 
     const TAxis* ax = fHist->GetXaxis();
     const TAxis* ay = fHist->GetYaxis();
     const double Nx = ax->GetNbins();
     const double Ny = ay->GetNbins();
 
     // Axes with limits where the user originally requested them, which we
     // adjusted to be the centres of the first and last bins.
     TH2* axes = new TH2C(UniqueName().c_str(),
                          TString::Format(";%s;%s",
                                          ax->GetTitle(), ay->GetTitle()),
                          Nx-1, BinCenter(ax, 1, fLogX), BinCenter(ax, Nx, fLogX),
                          Ny-1, BinCenter(ay, 1, fLogY), BinCenter(ay, Ny, fLogY));
     axes->Draw();
 
     if(fHist){
       // "colz same" will reuse axis's min and max, so set them helpfully here
       axes->SetMinimum(fHist->GetMinimum());
       axes->SetMaximum(fHist->GetMaximum());
     }
 
     axes->SetTitle(fHist->GetTitle());
     axes->GetXaxis()->SetLabelSize(ax->GetLabelSize());
     axes->GetYaxis()->SetLabelSize(ay->GetLabelSize());
     axes->GetXaxis()->SetLabelOffset(ax->GetLabelOffset());
     axes->GetYaxis()->SetLabelOffset(ay->GetLabelOffset());
     axes->GetXaxis()->CenterTitle();
     axes->GetYaxis()->CenterTitle();
 
     if(fLogX) gPad->SetLogx();
     if(fLogY) gPad->SetLogy();
 
     gPad->Update();
   }

void ana::Surface::FillSurface	(	const std::string &	progTitle,
		const IExperiment *	expt,
		osc::IOscCalcAdjustable *	calc,
		const FitAxis &	xax,
		const FitAxis &	yax,
		const std::vector< const IFitVar * > &	profVars,
		const std::vector< const ISyst * > &	profSysts,
		const std::map< const IFitVar *, std::vector< double >> &	seedPts,
		const std::vector< SystShifts > &	systSeedPts
	)

protected

Definition at line 145 of file sbnana/sbnana/CAFAna/Analysis/Surface.cxx.

   {
     if(fParallel && !(profVars.empty() && profSysts.empty())){
       // Minuit calls this at some point, and it creates a static. Which
       // doesn't like happening on multiple threads at once. So do it upfront
       // while we're still single-threaded.
       ROOT::Minuit2::StackAllocatorHolder::Get();
     }
 
     // Nothing created during surface filling belongs in a
     // directory. Unfortunately the local guards in Spectrum etc are racey when
     // run in parallel. But this should cover the whole lot safely.
     DontAddDirectory guard;
 
     Progress* prog = 0;
     // Difficult to keep a progress bar properly up to date when threaded
     if(!fParallel) prog = new Progress(progTitle);
     ThreadPool* pool = 0;
     if(fParallel){
       pool = new ThreadPool;
       pool->ShowProgress(progTitle);
     }
 
     const int Nx = fHist->GetNbinsX();
     const int Ny = fHist->GetNbinsY();
 
     // Fill bins in "random" order so that the progress bar is accurate
     // straight away instead of possibly being misled by whatever atypical
     // points we start with. This step is a prime which guarantees we get every
     // cell.
     int step = 7919;
     // Very unlikely (Nx or Ny is a multiple of step), but just to be safe.
     if((Nx*Ny)%step == 0) step = 1;
 
     int bin = 0;
     int neval = 0;
 
     // Allow the surface to be parallelised across multiple jobs by splitting
     // up the full surface into patches, and only running bins that fall inside
     // a certain patch
     int grid_stride = 1;
     if(RunningOnGrid() && NumJobs() > 1){
       grid_stride = NumJobs();
       bin = (step*JobNumber()) % (Nx*Ny);
     }
 
     do{
       const int x = bin%Nx+1;
       const int y = bin/Nx+1;
 
       const double xv = BinCenter(fHist->GetXaxis(), x, xax.islog);
       const double yv = BinCenter(fHist->GetYaxis(), y, yax.islog);
 
       // For parallel running need to set these at point of use otherwise we
       // race.
       if(!fParallel){
         xax.var->SetValue(calc, xv);
         yax.var->SetValue(calc, yv);
       }
 
       if(xax.var->Penalty(xv, calc) > 1e-10){
         std::cerr << "Warning! " << xax.var->ShortName() << " = " << xv
                   << " has penalty of " << xax.var->Penalty(xv, calc)
                   << " that could have been applied in surface. "
                   << "This should never happen." << std::endl;
       }
       if(yax.var->Penalty(yv, calc) > 1e-10){
         std::cerr << "Warning! " << yax.var->ShortName() << " = " << yv
                   << " has penalty of " << yax.var->Penalty(yv, calc)
                   << " that could have been applied in surface. "
                   << "This should never happen." << std::endl;
       }
 
       if(fParallel){
         pool->AddMemberTask(this, &Surface::FillSurfacePoint,
                             expt, calc,
                             xax, xv, yax, yv,
                             profVars, profSysts, seedPts, systSeedPts);
       }
       else{
         FillSurfacePoint(expt, calc,
                          xax, xv, yax, yv,
                          profVars, profSysts, seedPts, systSeedPts);
         ++neval;
         prog->SetProgress(neval*grid_stride/double(Nx*Ny));
       }
 
       if(grid_stride > 1) fBinMask.push_back(bin);
 
       for(int i = 0; i < grid_stride; ++i){
         bin = (bin + step) % (Nx * Ny);
         if(bin == 0) break;
       }
     } while(bin != 0);
 
 
     if(fParallel){
       pool->Finish();
       delete pool;
     }
     else{
       prog->Done();
       delete prog;
     }
   }

void ana::Surface::FillSurfacePoint	(	const IExperiment *	expt,
		osc::IOscCalcAdjustable *	calc,
		const FitAxis &	xax,
		double	x,
		const FitAxis &	yax,
		double	y,
		const std::vector< const IFitVar * > &	profVars,
		const std::vector< const ISyst * > &	profSysts,
		const std::map< const IFitVar *, std::vector< double >> &	seedPts,
		const std::vector< SystShifts > &	systSeedPts
	)

protected

Definition at line 291 of file sbnana/sbnana/CAFAna/Analysis/Surface.cxx.

   {
     osc::IOscCalcAdjustable* calcNoHash = 0; // specific to parallel mode
 
     if(fParallel){
       // Need to take our own copy so that we don't get overwritten by someone
       // else's changes.
       calc = calc->Copy();
       xax.var->SetValue(calc, x);
       yax.var->SetValue(calc, y);
 
       calcNoHash = new OscCalcNoHash(calc);
     }
 
     //Make sure that the profiled values of fitvars do not persist between steps.
     for(int i = 0; i < (int)fSeedValues.size(); ++i) profVars[i]->SetValue( calc, fSeedValues[i] );
 
     double chi;
     if(profVars.empty() && profSysts.empty()){
       chi = expt->ChiSq(fParallel ? calcNoHash : calc);
     }
     else{
       Fitter fitter(expt, profVars, profSysts);
       fitter.SetPrecision(fPrec);
       SystShifts bestSysts;
       chi = fitter.Fit(calc, bestSysts, seedPts, systSeedPts, Fitter::kQuiet);
 
       for(unsigned int i = 0; i < profVars.size(); ++i){
         fProfHists[i]->Fill(x, y, profVars[i]->GetValue(calc));
       }
       for(unsigned int j = 0; j < profSysts.size(); ++j){
         fProfHists[j+profVars.size()]->Fill(x, y, bestSysts.GetShift(profSysts[j]));
       }
     }
 
     fHist->Fill(x, y, chi);
 
     if(fParallel){
       delete calc;
       delete calcNoHash;
     }
   }

std::vector< TGraph * > ana::Surface::GetGraphs	(	TH2 *	fc,
		double	minchi = `-1`
	)

For expert use, custom painting of contours.

Definition at line 423 of file sbnana/sbnana/CAFAna/Analysis/Surface.cxx.

   {
     CheckMask("GetGraphs");
 
     std::vector<TGraph*> ret;
 
     if(minchi < 0) minchi = fMinChi;
     std::unique_ptr<TH2> surf((TH2*)fHist->Clone(UniqueName().c_str()));
     surf->Add(fc, -1);
 
     TVirtualPad* bak = gPad;
 
     const bool wasbatch = gROOT->IsBatch();
     gROOT->SetBatch(); // User doesn't want to see our temporary canvas
     TCanvas tmp;
 
     gStyle->SetOptStat(0);
 
     const double level = minchi-fMinChi;
     surf->SetContour(1, &level);
     surf->Draw("cont list");
 
     tmp.Update();
     tmp.Paint();
 
     gROOT->SetBatch(wasbatch);
     gPad = bak;
 
     // The graphs we need (contained inside TLists, contained inside
     // TObjArrays) are in the list of specials. But we need to be careful about
     // types, because other stuff can get in here too (TDatabasePDG for
     // example).
     TCollection* specs = gROOT->GetListOfSpecials();
 
     TIter nextSpec(specs);
     while(TObject* spec = nextSpec()){
       if(!spec->InheritsFrom(TObjArray::Class())) continue;
       TObjArray* conts = (TObjArray*)spec;
 
       if(conts->IsEmpty()) continue;
 
       if(!conts->At(0)->InheritsFrom(TList::Class())) continue;
       TList* cont = (TList*)conts->At(0);
 
       TIter nextObj(cont);
       // Contour could be split into multiple pieces
       while(TObject* obj = nextObj()){
         if(!obj->InheritsFrom(TGraph::Class())) continue;
 
         ret.push_back((TGraph*)obj->Clone(UniqueName().c_str()));
       } // end for obj
     } // end for spec
 
     return ret;
   }

std::vector<TH2*> ana::Surface::GetMarginalizedHists ( )

inline

Deprecated. Retained for backwards compatibility.

Definition at line 98 of file sbnana/sbnana/CAFAna/Analysis/Surface.h.

98 {return fProfHists;}

ana::Surface::fProfHists

std::vector< TH2 * > fProfHists

Definition: sbnana/sbnana/CAFAna/Analysis/Surface.h:139

double ana::Surface::GetMinX ( ) const

inline

Definition at line 82 of file sbnana/sbnana/CAFAna/Analysis/Surface.h.

82 {return fMinX;}

ana::Surface::fMinX

double fMinX

Definition: sbnana/sbnana/CAFAna/Analysis/Surface.h:136

double ana::Surface::GetMinY ( ) const

inline

Definition at line 83 of file sbnana/sbnana/CAFAna/Analysis/Surface.h.

83 {return fMinY;}

ana::Surface::fMinY

double fMinY

Definition: sbnana/sbnana/CAFAna/Analysis/Surface.h:136

std::vector<TH2*> ana::Surface::GetProfiledHists ( )

inline

Maps of the values taken on by the profiled parameters.

Definition at line 96 of file sbnana/sbnana/CAFAna/Analysis/Surface.h.

96 {return fProfHists;}

ana::Surface::fProfHists

std::vector< TH2 * > fProfHists

Definition: sbnana/sbnana/CAFAna/Analysis/Surface.h:139

std::unique_ptr< Surface > ana::Surface::LoadFrom ( TDirectory * dir )

static

Definition at line 597 of file sbnana/sbnana/CAFAna/Analysis/Surface.cxx.

   {
     DontAddDirectory guard;
 
     TObjString* tag = (TObjString*)dir->Get("type");
     assert(tag);
     assert(tag->GetString() == "Surface");
 
     const TVectorD v = *(TVectorD*)dir->Get("minValues");
     const TVectorD s = *(TVectorD*)dir->Get("seeds");
     const TVectorD* m = (TVectorD*)dir->Get("mask");
 
     std::unique_ptr<Surface> surf(new Surface);
     surf->fHist = (TH2F*)dir->Get("hist");
     surf->fMinChi = v[0];
     surf->fMinX = v[1];
     surf->fMinY = v[2];
 
     for(int idx = 0; idx < s.GetNrows(); ++idx){
       surf->fSeedValues.push_back(s[idx]);
     }
 
     for(int idx = 0; ; ++idx){
       TH2* h = (TH2*)dir->Get(TString::Format("profHists/hist%d", idx));
       if(h) surf->fProfHists.push_back(h); else break;
     }
 
     if(m){
       for(int idx = 0; idx < m->GetNrows(); ++idx)
         surf->fBinMask.push_back((*m)[idx]);
     }
 
     const TObjString* logx = (TObjString*)dir->Get("logx");
     const TObjString* logy = (TObjString*)dir->Get("logy");
     // Tolerate missing log tags for backwards compatibility
     surf->fLogX = (logx && logx->GetString() == "yes");
     surf->fLogY = (logy && logy->GetString() == "yes");
 
     return surf;
   }

std::unique_ptr< Surface > ana::Surface::LoadFromMulti	(	const std::vector< TFile * > &	files,
		const std::string &	label
	)

static

Definition at line 640 of file sbnana/sbnana/CAFAna/Analysis/Surface.cxx.

   {
     std::vector<std::unique_ptr<Surface>> surfs;
     for(TFile* f: files) {
       surfs.push_back(Surface::LoadFrom(f->GetDirectory(label.c_str())));
     }
 
     int Nx = surfs[0]->fHist->GetNbinsX();
     int Ny = surfs[0]->fHist->GetNbinsY();
     size_t nbins = Nx * Ny;
     std::vector<int> binMask;
 
     // Loop over the surfaces to check all bins are present
     for(size_t i = 0; i < surfs.size(); ++i) {
       for(int bin : surfs[i]->fBinMask) {
         if (std::count(binMask.begin(), binMask.end(), bin))
           assert(false && "Repeated bin found in surfaces being merged.");
         binMask.push_back(bin);
       }
     }
     if(binMask.size() != nbins) {
       std::cout << "Missing bins found in surfaces being merged. "
                 << "Are you sure you included all files for this surface?"
                 << std::endl;
       assert(false);
     }
 
     // Create return surface and initialise with first in list
     std::unique_ptr<Surface> ret(new Surface);
     ret->fHist = new TH2F(*(TH2F*)surfs[0]->ToTH2(0));
     ret->fMinChi = surfs[0]->fMinChi;
     ret->fMinX = surfs[0]->fMinX;
     ret->fMinY = surfs[0]->fMinY;
     ret->fLogX = surfs[0]->fLogX;
     ret->fLogY = surfs[0]->fLogY;
     ret->fSeedValues = surfs[0]->fSeedValues;
     for(TH2* h: surfs[0]->fProfHists)
       ret->fProfHists.push_back((TH2F*)h->Clone());
 
     // Loop over other surfaces and add them in
     for(size_t i = 1; i < surfs.size(); ++i) {
       ret->fHist->Add(surfs[i]->ToTH2(0));
       if(surfs[i]->fMinChi < ret->fMinChi) {
         ret->fMinChi = surfs[i]->fMinChi;
         ret->fMinX = surfs[i]->fMinX;
         ret->fMinY = surfs[i]->fMinY;
       }
       for(size_t j = 0; j < ret->fProfHists.size(); ++j)
         ret->fProfHists[j]->Add(surfs[i]->fProfHists[j]);
     }
 
     // Scale hist by global minimum
     for(int x = 0; x < ret->fHist->GetNbinsX()+2; ++x){
       for(int y = 0; y < ret->fHist->GetNbinsY()+2; ++y){
         ret->fHist->SetBinContent(x, y, ret->fHist->GetBinContent(x, y)-ret->fMinChi);
       }
     }
 
     ret->fHist->SetMinimum(0);
 
     return ret;
   }

std::unique_ptr< Surface > ana::Surface::LoadFromMulti	(	const std::string &	wildcard,
		const std::string &	label
	)

static

Definition at line 705 of file sbnana/sbnana/CAFAna/Analysis/Surface.cxx.

   {
     std::vector<std::string> fileNames = Wildcard(wildcard);
     std::vector<TFile*> files;
     for(std::string fileName : fileNames){
       files.push_back(TFile::Open(fileName.c_str()));
     }
     std::unique_ptr<Surface> ret = LoadFromMulti(files, label);
     for(TFile* f: files) delete f;
     return ret;
   }

double ana::Surface::MinChi ( ) const

inline

Definition at line 81 of file sbnana/sbnana/CAFAna/Analysis/Surface.h.

81 {return fMinChi;}

ana::Surface::fMinChi

double fMinChi

Definition: sbnana/sbnana/CAFAna/Analysis/Surface.h:135

void ana::Surface::SaveTo ( TDirectory * dir ) const

Definition at line 559 of file sbnana/sbnana/CAFAna/Analysis/Surface.cxx.

   {
     TDirectory* tmp = gDirectory;
     dir->cd();
     TObjString("Surface").Write("type");
 
     TVectorD v(3);
     v[0] = fMinChi;
     v[1] = fMinX;
     v[2] = fMinY;
     v.Write("minValues");
 
     fHist->Write("hist");
 
     TVectorD s(fSeedValues.size(), &fSeedValues[0]);
     s.Write("seeds");
 
     TDirectory* profDir = dir->mkdir("profHists");
     int idx = 0;
     for(auto it: fProfHists){
       profDir->cd();
       it->Write( TString::Format("hist%d", idx++));
     }
 
     if(!fBinMask.empty()){
       const std::vector<double> tmp(fBinMask.begin(), fBinMask.end());
       TVectorD m(tmp.size(), &tmp[0]);
       m.Write("mask");
     }
 
     dir->cd();
     TObjString(fLogX ? "yes" : "no").Write("logx");
     TObjString(fLogY ? "yes" : "no").Write("logy");
 
     tmp->cd();
   }

void ana::Surface::SetTitle ( const char * str )

Definition at line 518 of file sbnana/sbnana/CAFAna/Analysis/Surface.cxx.

   {
     fHist->SetTitle(str);
   }

TH2F * ana::Surface::ToTH2 ( double minchi = -1 ) const

Definition at line 499 of file sbnana/sbnana/CAFAna/Analysis/Surface.cxx.

   {
     // Could have a file temporarily open
     DontAddDirectory guard;
 
     TH2F* ret = new TH2F(*fHist);
 
     if(minchi >= 0){
       for(int x = 0; x < ret->GetNbinsX()+2; ++x){
         for(int y = 0; y < ret->GetNbinsY()+2; ++y){
           ret->SetBinContent(x, y, ret->GetBinContent(x, y)+fMinChi-minchi);
         }
       }
     }
 
     return ret;
   }