From af25ad18997a4c5a908b9f53936b2ebb727750f7 Mon Sep 17 00:00:00 2001 From: Andrew Gilbert Date: Mon, 17 Nov 2014 13:17:34 +0100 Subject: [PATCH 1/7] Cherry picking commit to SLC6 branch can clear the _morphParams pointers so these are updated to the variables copied into the workspace --- interface/VerticalInterpHistPdf.h | 1 + src/VerticalInterpHistPdf.cc | 6 ++++++ 2 files changed, 7 insertions(+) diff --git a/interface/VerticalInterpHistPdf.h b/interface/VerticalInterpHistPdf.h index ff16f696cc6..d2c518bf8da 100644 --- a/interface/VerticalInterpHistPdf.h +++ b/interface/VerticalInterpHistPdf.h @@ -285,6 +285,7 @@ class FastVerticalInterpHistPdf2Base : public RooAbsPdf { // initialize the morphParams and the sentry. to be called by the daughter class, sets also _initBase to true void initBase() const ; + virtual Bool_t importWorkspaceHook(RooWorkspace& ws); private: ClassDef(FastVerticalInterpHistPdf2Base,1) // diff --git a/src/VerticalInterpHistPdf.cc b/src/VerticalInterpHistPdf.cc index a1636609d06..7cc7458831b 100644 --- a/src/VerticalInterpHistPdf.cc +++ b/src/VerticalInterpHistPdf.cc @@ -747,6 +747,12 @@ FastVerticalInterpHistPdf2Base::FastVerticalInterpHistPdf2Base(const FastVertica } +Bool_t FastVerticalInterpHistPdf2Base::importWorkspaceHook(RooWorkspace& ws) { + _initBase = false; + _morphParams.clear(); + _sentry.reset(); + return kFALSE; +} //_____________________________________________________________________________ From 99e33a40955c7497842adef1aa930cd4d38cdac1 Mon Sep 17 00:00:00 2001 From: adavidzh Date: Thu, 20 Nov 2014 16:23:27 +0100 Subject: [PATCH 2/7] Bug fix for #138. Closes #153. --- interface/CloseCoutSentry.h | 2 ++ interface/utils.h | 5 ++-- src/CascadeMinimizer.cc | 9 ++++--- src/CloseCoutSentry.cc | 10 ++++++++ src/utils.cc | 51 +++++++++++++++++++++++++------------ 5 files changed, 55 insertions(+), 22 deletions(-) diff --git a/interface/CloseCoutSentry.h b/interface/CloseCoutSentry.h index f05a4c46c60..97538d71a01 100644 --- a/interface/CloseCoutSentry.h +++ b/interface/CloseCoutSentry.h @@ -13,6 +13,7 @@ class CloseCoutSentry { // break through any sentry, even the ones above myself (for critical error messages, or debug) static void breakFree() ; FILE *trueStdOut(); + static FILE *trueStdOutGlobal(); private: bool silent_; static int fdOut_, fdErr_; @@ -20,6 +21,7 @@ class CloseCoutSentry { // always clear, even if I was not the one closing it void static reallyClear() ; static FILE *trueStdOut_; + static CloseCoutSentry *owner_; bool stdOutIsMine_; }; diff --git a/interface/utils.h b/interface/utils.h index 651e5daa6ca..61c2eb71f6a 100644 --- a/interface/utils.h +++ b/interface/utils.h @@ -3,6 +3,7 @@ #include #include +#include struct RooDataHist; struct RooAbsData; struct RooAbsPdf; @@ -93,8 +94,8 @@ namespace utils { // Set range of physics model parameters void setModelParameterRanges( const std::string & setPhysicsModelParameterRangeExpression, const RooArgSet & params); - bool checkParameterBoundary( const RooRealVar &); - bool checkParameterBoundaries( const RooArgSet &); + bool isParameterAtBoundary( const RooRealVar &); + bool anyParameterAtBoundaries( const RooArgSet &, int verbosity); void reorderCombinations(std::vector > &, const std::vector &, const std::vector &); std::vector > generateCombinations(const std::vector &vec); diff --git a/src/CascadeMinimizer.cc b/src/CascadeMinimizer.cc index b71a85899de..da7b5fcf5b6 100644 --- a/src/CascadeMinimizer.cc +++ b/src/CascadeMinimizer.cc @@ -278,10 +278,11 @@ bool CascadeMinimizer::minimize(int verbose, bool cascade) nllParams->remove(CascadeMinimizerGlobalConfigs::O().pdfCategories); nllParams->remove(CascadeMinimizerGlobalConfigs::O().parametersOfInterest); - bool boundariesOk = utils::checkParameterBoundaries(*nllParams); - if(!boundariesOk){ - std::cout << " [WARNING] After the fit some parameters are at their boundary (see above)." << std::endl; - std::cout << " [WARNING] Are you sure your model is correct?" << std::endl; + bool boundariesNotOk = utils::anyParameterAtBoundaries(*nllParams, verbose); + if(boundariesNotOk && verbose >= 1){ + fprintf(CloseCoutSentry::trueStdOutGlobal(), + " [WARNING] After the fit some parameters are at their boundary.\n" + " [WARNING] Are you sure your model is correct?\n"); } return ret; diff --git a/src/CloseCoutSentry.cc b/src/CloseCoutSentry.cc index ab239d5807b..47c545a4e70 100644 --- a/src/CloseCoutSentry.cc +++ b/src/CloseCoutSentry.cc @@ -7,6 +7,8 @@ bool CloseCoutSentry::open_ = true; int CloseCoutSentry::fdOut_ = 0; int CloseCoutSentry::fdErr_ = 0; FILE * CloseCoutSentry::trueStdOut_ = 0; +CloseCoutSentry *CloseCoutSentry::owner_ = 0; + CloseCoutSentry::CloseCoutSentry(bool silent) : silent_(silent), stdOutIsMine_(false) @@ -20,6 +22,7 @@ CloseCoutSentry::CloseCoutSentry(bool silent) : } freopen("/dev/null", "w", stdout); freopen("/dev/null", "w", stderr); + owner_ = this; } else { silent_ = false; } @@ -47,6 +50,7 @@ void CloseCoutSentry::reallyClear() sprintf(buf, "/dev/fd/%d", fdOut_); freopen(buf, "w", stdout); sprintf(buf, "/dev/fd/%d", fdErr_); freopen(buf, "w", stderr); open_ = true; + owner_ = 0; fdOut_ = fdErr_ = 0; } } @@ -56,6 +60,12 @@ void CloseCoutSentry::breakFree() reallyClear(); } +FILE *CloseCoutSentry::trueStdOutGlobal() +{ + if (!owner_) return stdout; + return owner_->trueStdOut(); +} + FILE *CloseCoutSentry::trueStdOut() { if (open_) return stdout; diff --git a/src/utils.cc b/src/utils.cc index d32209634bb..8896edc292f 100644 --- a/src/utils.cc +++ b/src/utils.cc @@ -6,6 +6,7 @@ #include #include #include +#include #include #include #include @@ -414,9 +415,11 @@ void utils::copyAttributes(const RooAbsArg &from, RooAbsArg &to) { if (!attribs.empty()) { for (std::set::const_iterator it = attribs.begin(), ed = attribs.end(); it != ed; ++it) to.setAttribute(it->c_str()); } - const std::map strattribs = from.stringAttributes(); + const std::map + strattribs = from.stringAttributes(); if (!strattribs.empty()) { - for (std::map::const_iterator it = strattribs.begin(), ed = strattribs.end(); it != ed; ++it) to.setStringAttribute(it->first.c_str(), it->second.c_str()); + for (std::map + ::const_iterator it = strattribs.begin(), ed = strattribs.end(); it != ed; ++it) to.setStringAttribute(it->first.c_str(), it->second.c_str()); } } @@ -695,7 +698,8 @@ std::vector > utils::generateCombinations(const std::vector timesFoundAtBoundary; + bool isAnyBad = true; RooLinkedListIter iter = params.iterator(); int i = 0; for (RooRealVar *a = (RooRealVar *) iter.Next(); a != 0; a = (RooRealVar *) iter.Next(), ++i) { - bool isBad = checkParameterBoundary(*a); - isNoneBad &= isBad; + + bool isBad = isParameterAtBoundary(*a); + + if(isBad){ + std::string varName((*a).GetName()); + + if( verbosity >= 9 || (timesFoundAtBoundary[varName] < 3 && verbosity > -1) ){ + fprintf(CloseCoutSentry::trueStdOutGlobal()," [WARNING] Found [%s] at boundary. \n", (*a).GetName()); + std::cout << " "; (*a).Print(); + } + + timesFoundAtBoundary[varName]++; + } + + isAnyBad |= isBad; } - return isNoneBad; + // for( std::unordered_map::value_type e : timesFoundAtBoundary ){ + // printf("e %s -> %i\n", e.first.c_str(), e.second); + // } + + return isAnyBad; } - From bdac319bbb54bb5b8460e24e5c72079d6008610f Mon Sep 17 00:00:00 2001 From: nick Date: Thu, 20 Nov 2014 15:58:38 +0100 Subject: [PATCH 3/7] Adding plots to output of diffNuisances if output file specified, additional plots of the post fit and prefit nuisances and errors are saved. --- test/diffNuisances.py | 88 ++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 87 insertions(+), 1 deletion(-) diff --git a/test/diffNuisances.py b/test/diffNuisances.py index 0783014d19e..ac4c710ed1f 100644 --- a/test/diffNuisances.py +++ b/test/diffNuisances.py @@ -46,6 +46,12 @@ fpf_b = fit_b.floatParsFinal() fpf_s = fit_s.floatParsFinal() pulls = [] + +nuis_p_i=0 +# Also make histograms for pull distributions: +hist_fit_b = ROOT.TH1F("prefit_fit_b" ,"B-only fit Nuisances;;#theta ",prefit.getSize(),0,prefit.getSize()) +hist_fit_s = ROOT.TH1F("prefit_fit_s" ,"S+B fit Nuisances ;;#theta ",prefit.getSize(),0,prefit.getSize()) +hist_prefit = ROOT.TH1F("prefit_nuisancs","Prefit Nuisances ;;#theta ",prefit.getSize(),0,prefit.getSize()) for i in range(fpf_s.getSize()): nuis_s = fpf_s.at(i) name = nuis_s.GetName(); @@ -65,7 +71,22 @@ row += [ " n/a " ] else: row += [ "%+.2f +/- %.2f" % (nuis_x.getVal(), nuis_x.getError()) ] + if nuis_p != None: + if options.plotfile: + if fit_name=='b': + nuis_p_i+=1 + hist_fit_b.SetBinContent(nuis_p_i,nuis_x.getVal()) + hist_fit_b.SetBinError(nuis_p_i,nuis_x.getError()) + hist_fit_b.GetXaxis().SetBinLabel(nuis_p_i,name) + if fit_name=='s': + hist_fit_s.SetBinContent(nuis_p_i,nuis_x.getVal()) + hist_fit_s.SetBinError(nuis_p_i,nuis_x.getError()) + hist_fit_s.GetXaxis().SetBinLabel(nuis_p_i,name) + hist_prefit.SetBinContent(nuis_p_i,mean_p) + hist_prefit.SetBinError(nuis_p_i,sigma_p) + hist_prefit.GetXaxis().SetBinLabel(nuis_p_i,name) + valShift = (nuis_x.getVal() - mean_p)/sigma_p if fit_name == 'b': pulls.append(valShift) @@ -167,6 +188,7 @@ if options.plotfile: import ROOT + fout = ROOT.TFile(options.plotfile,"RECREATE") ROOT.gROOT.SetStyle("Plain") ROOT.gStyle.SetOptFit(1) histogram = ROOT.TH1F("pulls", "Pulls", 60, -3, 3) @@ -179,4 +201,68 @@ histogram.SetMarkerSize(2) #histogram.Fit("gaus") histogram.Draw("pe") - canvas.SaveAs(options.plotfile) + #canvas.SaveAs(options.plotfile) + fout.WriteTObject(canvas) + + canvas_nuis = ROOT.TCanvas("nuisancs", "nuisances", 900, 600) + hist_fit_s.SetLineColor(ROOT.kRed) + hist_fit_s.SetMarkerColor(ROOT.kRed) + hist_fit_b.SetLineColor(ROOT.kBlue) + hist_fit_b.SetMarkerColor(ROOT.kBlue) + hist_fit_b.SetMarkerStyle(20) + hist_fit_s.SetMarkerStyle(20) + hist_fit_b.SetMarkerSize(1.0) + hist_fit_s.SetMarkerSize(1.0) + hist_fit_b.SetLineWidth(2) + hist_fit_s.SetLineWidth(2) + hist_prefit.SetLineWidth(2) + hist_prefit.SetTitle("Nuisance Paramaeters") + hist_prefit.SetLineColor(ROOT.kBlack) + hist_prefit.SetFillColor(ROOT.kGray) + hist_prefit.Draw("E2") + hist_prefit.Draw("histsame") + hist_fit_b.Draw("E1Psame") + hist_fit_s.Draw("E1Psame") + canvas_nuis.RedrawAxis() + leg=ROOT.TLegend(0.6,0.7,0.89,0.89) + leg.SetFillColor(0) + leg.SetTextFont(42) + leg.AddEntry(hist_prefit,"Prefit","FL") + leg.AddEntry(hist_fit_b,"B-only fit","EPL") + leg.AddEntry(hist_fit_s,"S+B fit" ,"EPL") + leg.Draw() + fout.WriteTObject(canvas_nuis) + canvas_pferrs = ROOT.TCanvas("post_fit_errs", "post_fit_errs", 900, 600) + hist_fit_e_s = hist_fit_s.Clone() + hist_fit_e_b = hist_fit_b.Clone() + for b in range(1,hist_fit_e_s.GetNbinsX()+1): + hist_fit_e_s.SetBinContent(b,hist_fit_s.GetBinError(b)/hist_prefit.GetBinError(b)) + hist_fit_e_b.SetBinContent(b,hist_fit_b.GetBinError(b)/hist_prefit.GetBinError(b)) + hist_fit_e_s.SetBinError(b,0) + hist_fit_e_b.SetBinError(b,0) + hist_fit_e_s.SetFillColor(ROOT.kRed) + hist_fit_e_b.SetFillColor(ROOT.kBlue) + hist_fit_e_s.SetBarWidth(0.4) + hist_fit_e_b.SetBarWidth(0.4) + hist_fit_e_b.SetBarOffset(0.45) + hist_fit_e_b.GetYaxis().SetTitle("#sigma_{#theta}/(#sigma_{#theta} prefit)") + hist_fit_e_b.SetTitle("Nuisance Parameter Uncertainty Reduction") + hist_fit_e_b.SetMaximum(1.5) + hist_fit_e_b.SetMinimum(0) + hist_fit_e_b.Draw("bar") + hist_fit_e_s.Draw("barsame") + leg_rat=ROOT.TLegend(0.6,0.7,0.89,0.89) + leg_rat.SetFillColor(0) + leg_rat.SetTextFont(42) + leg_rat.AddEntry(hist_fit_e_b,"B-only fit","F") + leg_rat.AddEntry(hist_fit_e_s,"S+B fit" ,"F") + leg_rat.Draw() + line_one = ROOT.TLine(0,1,hist_fit_e_s.GetXaxis().GetXmax(),1) + line_one.SetLineColor(1); line_one.SetLineStyle(2); line_one.SetLineWidth(2) + line_one.Draw() + canvas_pferrs.RedrawAxis() + + fout.WriteTObject(canvas_pferrs) + + + From 5062fa53561c95daf26336e2bbda3469a69db1ce Mon Sep 17 00:00:00 2001 From: adavidzh Date: Thu, 20 Nov 2014 17:36:15 +0100 Subject: [PATCH 4/7] Invert initaliser after having inverted logic. --- src/utils.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils.cc b/src/utils.cc index 8896edc292f..c5db9ad5365 100644 --- a/src/utils.cc +++ b/src/utils.cc @@ -722,7 +722,7 @@ bool utils::isParameterAtBoundary( const RooRealVar ¶m ){ bool utils::anyParameterAtBoundaries( const RooArgSet ¶ms, int verbosity ){ static std::unordered_map timesFoundAtBoundary; - bool isAnyBad = true; + bool isAnyBad = false; RooLinkedListIter iter = params.iterator(); int i = 0; for (RooRealVar *a = (RooRealVar *) iter.Next(); a != 0; a = (RooRealVar *) iter.Next(), ++i) { From 60a62a2b1b702060b4577db9286b21cb7e7d0135 Mon Sep 17 00:00:00 2001 From: nucleosynthesis Date: Sat, 22 Nov 2014 00:57:07 +0100 Subject: [PATCH 5/7] fix post-fit nuisance plots --- test/diffNuisances.py | 31 +++++++++++++++++++++++++------ 1 file changed, 25 insertions(+), 6 deletions(-) diff --git a/test/diffNuisances.py b/test/diffNuisances.py index ac4c710ed1f..31a2504059f 100644 --- a/test/diffNuisances.py +++ b/test/diffNuisances.py @@ -186,6 +186,18 @@ elif options.format == "html": print "" +def getGraph(hist,shift): + + gr = ROOT.TGraphErrors() + gr.SetName(hist.GetName()) + for i in range(hist.GetNbinsX()): + x = hist.GetBinCenter(i+1)+shift + y = hist.GetBinContent(i+1) + e = hist.GetBinError(i+1) + gr.SetPoint(i,x,y) + gr.SetPointError(i,float(abs(shift))*0.8,e) + return gr + if options.plotfile: import ROOT fout = ROOT.TFile(options.plotfile,"RECREATE") @@ -205,6 +217,10 @@ fout.WriteTObject(canvas) canvas_nuis = ROOT.TCanvas("nuisancs", "nuisances", 900, 600) + hist_fit_e_s = hist_fit_s.Clone() + hist_fit_e_b = hist_fit_b.Clone() + hist_fit_s = getGraph(hist_fit_s,-0.1) + hist_fit_b = getGraph(hist_fit_b, 0.1) hist_fit_s.SetLineColor(ROOT.kRed) hist_fit_s.SetMarkerColor(ROOT.kRed) hist_fit_b.SetLineColor(ROOT.kBlue) @@ -219,11 +235,15 @@ hist_prefit.SetTitle("Nuisance Paramaeters") hist_prefit.SetLineColor(ROOT.kBlack) hist_prefit.SetFillColor(ROOT.kGray) + hist_prefit.SetMaximum(3) + hist_prefit.SetMinimum(-3) hist_prefit.Draw("E2") hist_prefit.Draw("histsame") - hist_fit_b.Draw("E1Psame") - hist_fit_s.Draw("E1Psame") + hist_fit_b.Draw("EPsame") + hist_fit_s.Draw("EPsame") + canvas_nuis.SetGridx() canvas_nuis.RedrawAxis() + canvas_nuis.RedrawAxis('g') leg=ROOT.TLegend(0.6,0.7,0.89,0.89) leg.SetFillColor(0) leg.SetTextFont(42) @@ -232,12 +252,11 @@ leg.AddEntry(hist_fit_s,"S+B fit" ,"EPL") leg.Draw() fout.WriteTObject(canvas_nuis) + canvas_pferrs = ROOT.TCanvas("post_fit_errs", "post_fit_errs", 900, 600) - hist_fit_e_s = hist_fit_s.Clone() - hist_fit_e_b = hist_fit_b.Clone() for b in range(1,hist_fit_e_s.GetNbinsX()+1): - hist_fit_e_s.SetBinContent(b,hist_fit_s.GetBinError(b)/hist_prefit.GetBinError(b)) - hist_fit_e_b.SetBinContent(b,hist_fit_b.GetBinError(b)/hist_prefit.GetBinError(b)) + hist_fit_e_s.SetBinContent(b,hist_fit_e_s.GetBinError(b)/hist_prefit.GetBinError(b)) + hist_fit_e_b.SetBinContent(b,hist_fit_e_b.GetBinError(b)/hist_prefit.GetBinError(b)) hist_fit_e_s.SetBinError(b,0) hist_fit_e_b.SetBinError(b,0) hist_fit_e_s.SetFillColor(ROOT.kRed) From bf6175e4b2c29a39f22eaecaa4d83420c580297a Mon Sep 17 00:00:00 2001 From: nick Date: Thu, 13 Nov 2014 17:08:43 +0100 Subject: [PATCH 6/7] Adding ability to group nuisances NOT in a specfic group run --freezeNuisanceGroup ^name to freeze all nuisances not in group "name" Additionally fixed a bug causing a crash when using freezeNuisanceGroup accessing the mc_bonly Some print out for when an unkwon group is asked for Conflicts: src/Combine.cc --- src/Combine.cc | 24 ++++++++++++++++++++++-- 1 file changed, 22 insertions(+), 2 deletions(-) diff --git a/src/Combine.cc b/src/Combine.cc index e33e5b9933d..897d4fec82b 100644 --- a/src/Combine.cc +++ b/src/Combine.cc @@ -449,14 +449,34 @@ void Combine::run(TString hlfFile, const std::string &dataset, double &limit, do std::vector nuisanceGroups; boost::algorithm::split(nuisanceGroups,freezeNuisanceGroups_,boost::algorithm::is_any_of(",")); for (std::vector::iterator ng_it=nuisanceGroups.begin();ng_it!=nuisanceGroups.end();ng_it++){ - RooArgSet toFreeze(*(w->set(Form("group_%s",(*ng_it).c_str())))); + bool freeze_complement=false; + if (boost::algorithm::starts_with((*ng_it),"^")){ + freeze_complement=true; + (*ng_it).erase(0,1); + } + + if (! w->set(Form("group_%s",(*ng_it).c_str()))){ + std::cerr << "Unknown nuisance group: " << (*ng_it) << std::endl; + throw std::invalid_argument("Unknown nuisance group name"); + } + RooArgSet groupNuisances(*(w->set(Form("group_%s",(*ng_it).c_str())))); + RooArgSet toFreeze; + + if (freeze_complement) { + RooArgSet still_floating(*mc->GetNuisanceParameters()); + still_floating.remove(groupNuisances,true,true); + toFreeze.add(still_floating); + } else { + toFreeze.add(groupNuisances); + } + if (verbose > 0) std::cout << "Freezing the following nuisance parameters: "; toFreeze.Print(""); utils::setAllConstant(toFreeze, true); if (nuisances) { RooArgSet newnuis(*nuisances); newnuis.remove(toFreeze, /*silent=*/true, /*byname=*/true); mc->SetNuisanceParameters(newnuis); - mc_bonly->SetNuisanceParameters(newnuis); + if (mc_bonly) mc_bonly->SetNuisanceParameters(newnuis); nuisances = mc->GetNuisanceParameters(); } } From 53fcf82b2efc5639de40019c691a1b90100d8500 Mon Sep 17 00:00:00 2001 From: nick Date: Tue, 25 Nov 2014 13:48:48 +0100 Subject: [PATCH 7/7] no starts_with in boost::algorithm (gcc version) --- src/Combine.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Combine.cc b/src/Combine.cc index 897d4fec82b..862eaa9bc5c 100644 --- a/src/Combine.cc +++ b/src/Combine.cc @@ -450,7 +450,7 @@ void Combine::run(TString hlfFile, const std::string &dataset, double &limit, do boost::algorithm::split(nuisanceGroups,freezeNuisanceGroups_,boost::algorithm::is_any_of(",")); for (std::vector::iterator ng_it=nuisanceGroups.begin();ng_it!=nuisanceGroups.end();ng_it++){ bool freeze_complement=false; - if (boost::algorithm::starts_with((*ng_it),"^")){ + if (((*ng_it).at(0)==std::string("^"))){ freeze_complement=true; (*ng_it).erase(0,1); }