15 #include "TObjString.h" 18 #include "CAFAna/Core/Var.h" 36 const std::vector<const IFitVar*>& profVars,
37 const std::vector<const ISyst*>& profSysts,
39 const std::vector<SystShifts>& systSeedPts,
42 : fParallel(parallel), fFitOpts(opts)
51 FillSurface(expt, calc, xax.
var, yax.
var, profVars, profSysts, seedPts, systSeedPts);
53 FindMinimum(expt, calc, xax.
var, yax.
var, profVars, profSysts, seedPts, systSeedPts);
64 const std::vector<const IFitVar*>& profVars,
65 const std::vector<const ISyst*>& profSysts)
71 for(
unsigned int i = 0;
i < profVars.size()+profSysts.size(); ++
i){
73 if(
i < profVars.size())
74 title = profVars[
i]->LatexName();
76 title = profSysts[
i-profVars.size()]->LatexName();
88 const std::vector<const IFitVar*>& profVars,
89 const std::vector<const ISyst*>& profSysts)
const 93 if(!profVars.empty() || !profSysts.empty()){
94 title +=
" (profiling ";
96 for(
const IFitVar*
v: profVars) title +=
v->ShortName() +
", ";
97 for(
const ISyst*
s: profSysts) title +=
s->ShortName() +
", ";
100 title.resize(title.size()-2);
111 const std::vector<const IFitVar *> &profVars,
112 const std::vector<const ISyst *> &profSysts,
114 const std::vector<SystShifts> &systSeedPts)
140 const int Nx =
fHist->GetNbinsX();
141 const int Ny =
fHist->GetNbinsY();
149 if ((Nx * Ny) % step == 0) step = 1;
157 int first = 0, last = Nx * Ny;
162 last = (job+1) * stride;
168 bin = (bin +
step) % (Nx * Ny);
170 }
else if (neval >= last) {
174 const int x = bin % Nx + 1;
175 const int y = bin / Nx + 1;
180 if (xvar->
Penalty(xv, calc) > 1
e-10)
183 <<
" has penalty of " << xvar->
Penalty(xv, calc)
184 <<
" that could have been applied in surface. " 185 <<
"This should never happen." <<
std::endl;
187 if (yvar->
Penalty(yv, calc) > 1
e-10)
190 <<
" has penalty of " << yvar->
Penalty(yv, calc)
191 <<
" that could have been applied in surface. " 192 <<
"This should never happen." <<
std::endl;
198 profVars, profSysts, seedPts, systSeedPts);
207 prog->
SetProgress(
double(neval-first) /
double(last-first));
212 bin = (bin +
step) % (Nx * Ny);
231 const std::vector<const IFitVar*>& profVars,
232 const std::vector<const ISyst*>& profSysts,
234 const std::vector<SystShifts>& systSeedPts)
249 if(profVars.empty() && profSysts.empty()){
250 chi = expt->
ChiSq(calc);
258 for(
unsigned int i = 0; i < profVars.size(); ++
i){
259 fProfHists[
i]->Fill(x, y, profVars[i]->GetValue(calc));
261 for(
unsigned int j = 0;
j < profSysts.size(); ++
j){
266 fHist->Fill(x, y, chi);
278 const std::vector<const IFitVar*>& profVars,
279 const std::vector<const ISyst*>& profSysts,
281 const std::vector<SystShifts>& systSeedPts)
284 double minchi = 1e10;
285 int minx =
fHist->GetNbinsX()/2;
286 int miny =
fHist->GetNbinsY()/2;
287 for(
int x = 1;
x <=
fHist->GetNbinsX(); ++
x){
288 for(
int y = 1;
y <=
fHist->GetNbinsY(); ++
y){
290 const double chi =
fHist->GetBinContent(
x,
y);
301 std::vector<const IFitVar*> allVars = {xvar, yvar};
302 allVars.insert(allVars.end(), profVars.begin(), profVars.end());
314 for(
int x = 0;
x <
fHist->GetNbinsX()+2; ++
x){
315 for(
int y = 0;
y <
fHist->GetNbinsY()+2; ++
y){
320 fHist->SetMinimum(0);
326 TDirectory *
tmp = gDirectory;
328 dir = dir->mkdir(name.c_str());
331 TObjString(
"FrequentistSurface").Write(
"type");
335 TDirectory *profDir = dir->mkdir(
"profHists");
352 dir = dir->GetDirectory(name.c_str());
357 TObjString *
tag = (TObjString *) dir->Get(
"type");
359 assert(tag->GetString() ==
"FrequentistSurface" || tag->GetString() ==
"Surface");
365 for (std::size_t
idx = 0;
idx < surf->fSeedValues.size(); ++
idx)
370 surf->fProfHists.push_back(h);
372 surf->fProfHists.push_back((TH2 *) dir->Get(
TString::Format(
"margHists/hist%lu",
idx)));
375 const TObjString*
logx = (TObjString*)dir->Get(
"logx");
376 const TObjString*
logy = (TObjString*)dir->Get(
"logy");
378 surf->fLogX = (logx && logx->GetString() ==
"yes");
379 surf->fLogY = (logy && logy->GetString() ==
"yes");
390 std::vector<std::unique_ptr<FrequentistSurface>> surfs;
391 for (TFile*
f : files) {
395 int Nx = surfs[0]->fHist->GetNbinsX();
396 int Ny = surfs[0]->fHist->GetNbinsY();
397 size_t nbins = Nx * Ny;
398 std::vector<int> binMask;
401 for (
size_t i = 0;
i < surfs.size(); ++
i) {
404 assert(
false &&
"Repeated bin found in surfaces being merged.");
405 binMask.push_back(
bin);
408 if (binMask.size() !=
nbins) {
409 std::cout <<
"Missing bins found in surfaces being merged. " 410 <<
"Are you sure you included all files for this surface?" 419 ret->fHist =
new TH2F(*(TH2F*)surfs[0]->
ToTH2(0));
420 ret->fBestLikelihood = surfs[0]->fBestLikelihood;
421 ret->fBestFitX = surfs[0]->fBestFitX;
422 ret->fBestFitY = surfs[0]->fBestFitY;
423 ret->fLogX = surfs[0]->fLogX;
424 ret->fLogY = surfs[0]->fLogY;
425 ret->fSeedValues = surfs[0]->fSeedValues;
427 ret->fProfHists.push_back((TH2F*)
h->Clone());
430 for (
size_t i = 1;
i < surfs.size(); ++
i) {
431 ret->fHist->Add(surfs[
i]->
ToTH2(0));
433 ret->fBestLikelihood = surfs[
i]->fBestLikelihood;
434 ret->fBestFitX = surfs[
i]->fBestFitX;
435 ret->fBestFitY = surfs[
i]->fBestFitY;
437 for (
size_t j = 0;
j < ret->fProfHists.size(); ++
j)
438 ret->fProfHists[
j]->Add(surfs[
i]->fProfHists[
j]);
442 for(
int x = 0;
x < ret->fHist->GetNbinsX()+2; ++
x){
443 for(
int y = 0;
y < ret->fHist->GetNbinsY()+2; ++
y){
444 ret->fHist->SetBinContent(
x,
y, ret->fHist->GetBinContent(
x,
y)-ret->fBestLikelihood);
448 ret->fHist->SetMinimum(0);
457 std::vector<std::string> fileNames =
Wildcard(wildcard);
458 std::vector<TFile*>
files;
460 files.push_back(TFile::Open(
fileName.c_str(),
"read"));
463 for(TFile*
f: files)
delete f;
void Finish()
Wait for all threads to complete before returning.
const std::string & LatexName() const
virtual double ChiSq(osc::IOscCalcAdjustable *osc, const SystShifts &syst=SystShifts::Nominal()) const
TH2 * Gaussian90Percent1D(const FrequentistSurface &s)
Up-value surface for 90% confidence in 1D in gaussian approximation.
size_t JobNumber()
What's the process number for a grid job?
virtual _IOscCalcAdjustable< T > * Copy() const =0
Cuts and Vars for the 2020 FD DiF Study.
MinuitFitter::FitOpts fFitOpts
double BinCenterX(int bin) const
double FillSurfacePoint(const IExperiment *expt, osc::IOscCalcAdjustable *calc, const IFitVar *xvar, double x, const IFitVar *yvar, double y, const std::vector< const IFitVar * > &profVars, const std::vector< const ISyst * > &profSysts, const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts)
bool RunningOnGrid()
Is this a grid (condor) job?
std::function< void(void)> func_t
The type of the user's worker functions.
Simple record of shifts applied to systematic parameters.
TH2 * Gaussian68Percent2D(const FrequentistSurface &s)
Up-value surface for 68% confidence in 2D in gaussian approximation.
double BinCenterY(int bin) const
TH2 * Gaussian2Sigma2D(const FrequentistSurface &s)
Up-value surface for 2 sigma confidence in 2D in gaussian approximation.
TH2 * Gaussian2Sigma1D(const FrequentistSurface &s)
Up-value surface for 2 sigma confidence in 1D in gaussian approximation.
static SystShifts Nominal()
std::string ProgressBarTitle(const IFitVar *xvar, const IFitVar *yvar, const std::vector< const IFitVar * > &profVars, const std::vector< const ISyst * > &profSysts) const
Encapsulate code to systematically shift a caf::SRProxy.
TH2F * ExpandedHistogram(const std::string &title, int nbinsx, double xmin, double xmax, bool xlog, int nbinsy, double ymin, double ymax, bool ylog)
Internal helper for Surface and FCSurface.
std::vector< TH2 * > fProfHists
Log-likelihood scan across two parameters.
TH2 * Gaussian95Percent2D(const FrequentistSurface &s)
Up-value surface for 95% confidence in 2D in gaussian approximation.
std::vector< std::string > Wildcard(const std::string &wildcardString)
Find files matching a UNIX glob, plus expand environment variables.
void FindMinimum(const IExperiment *expt, osc::IOscCalcAdjustable *calc, const IFitVar *xvar, const IFitVar *yvar, const std::vector< const IFitVar * > &profVars, const std::vector< const ISyst * > &profSysts, const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts)
std::vector< std::string > wildcard(const std::string &wildcardString)
virtual std::unique_ptr< IFitSummary > Fit(osc::IOscCalcAdjustable *seed, SystShifts &bestSysts=junkShifts, const SeedList &seedPts=SeedList(), const std::vector< SystShifts > &systSeedPts={}, Verbosity verb=kVerbose) const
Master fitting method. Depends on FitHelper and FitHelperSeeded.
virtual double Penalty(double, osc::IOscCalcAdjustable *) const
TH2 * Gaussian3Sigma2D(const FrequentistSurface &s)
Up-value surface for 3 sigma confidence in 2D in gaussian approximation.
T GetShift(const ISyst *syst) const
TH2 * Gaussian99Percent2D(const FrequentistSurface &s)
Up-value surface for 99% confidence in 2D in gaussian approximation.
static std::unique_ptr< FrequentistSurface > LoadFromMulti(std::vector< TFile * > files, std::string label)
TH2 * Gaussian3Sigma1D(const FrequentistSurface &s)
Up-value surface for 3 sigma confidence in 1D in gaussian approximation.
void SaveToHelper(TDirectory *dir) const
dir should already be the appropriate sub-directory
void AddTask(const func_t &func)
std::vector< double > fSeedValues
TH2 * ToTH2(double minchi=-1) const
virtual ~FrequentistSurface()
void SetProgress(double frac)
Update the progress fraction between zero and one.
const std::string & ShortName() const
virtual void SetValue(osc::IOscCalcAdjustable *osc, double val) const =0
TH2 * Gaussian90Percent2D(const FrequentistSurface &s)
Up-value surface for 90% confidence in 2D in gaussian approximation.
void SaveTo(TDirectory *dir, const std::string &name) const
Base class defining interface for experiments.
friend TH2 * Flat(double, const ISurface &)
Helper function for the gaussian approximation surfaces.
std::vector< int > fBinMask
A very simple thread pool for use by Surface.
assert(nhit_max >=nhit_nbins)
A simple ascii-art progress bar.
void ShowProgress(const std::string &title)
Interface definition for fittable variables.
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
TH2 * Gaussian99Percent1D(const FrequentistSurface &s)
Up-value surface for 99% confidence in 1D in gaussian approximation.
static void FillSurfObj(ISurface &surf, TDirectory *dir)
void CreateHistograms(const FitAxis &xax, const FitAxis &yax, const std::vector< const IFitVar * > &profVars, const std::vector< const ISyst * > &profSysts)
Prevent histograms being added to the current directory.
TH2 * Gaussian95Percent1D(const FrequentistSurface &s)
Up-value surface for 95% confidence in 1D in gaussian approximation.
void Done()
Call this when action is completed.
static std::unique_ptr< FrequentistSurface > LoadFrom(TDirectory *dir, const std::string &name)
TH2 * Gaussian68Percent1D(const FrequentistSurface &s)
Up-value surface for 68% confidence in 1D in gaussian approximation.
Collect information describing the axis of a fit variable.
void SetFitOpts(FitOpts opts)
virtual double GetValue(const osc::IOscCalcAdjustable *osc) const =0
virtual void FillSurface(const IExperiment *expt, osc::IOscCalcAdjustable *calc, const IFitVar *xvar, const IFitVar *yvar, const std::vector< const IFitVar * > &profVars, const std::vector< const ISyst * > &profSysts, const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts)
Perform MINUIT fits in one or two dimensions.