CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Public Member Functions | Protected Member Functions | Protected Attributes | Static Protected Attributes
MaxLikelihoodFit Class Reference

#include <MaxLikelihoodFit.h>

Inheritance diagram for MaxLikelihoodFit:
FitterAlgoBase LimitAlgo

Public Member Functions

virtual void applyOptions (const boost::program_options::variables_map &vm)
 
 MaxLikelihoodFit ()
 
virtual const std::string & name () const
 
virtual void setNToys (const int)
 
virtual void setToyNumber (const int)
 
 ~MaxLikelihoodFit ()
 
- Public Member Functions inherited from FitterAlgoBase
void applyOptionsBase (const boost::program_options::variables_map &vm)
 
 FitterAlgoBase (const char *title="<FillMe> specific options")
 
virtual bool run (RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, double &limit, double &limitErr, const double *hint)
 
- Public Member Functions inherited from LimitAlgo
virtual void applyDefaultOptions ()
 
 LimitAlgo ()
 
 LimitAlgo (const char *desc)
 
const
boost::program_options::options_description & 
options () const
 
virtual ~LimitAlgo ()
 

Protected Member Functions

void createFitResultTrees (const RooStats::ModelConfig &)
 
void getNormalizations (RooAbsPdf *pdf, const RooArgSet &obs, RooArgSet &out)
 
virtual bool runSpecific (RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, double &limit, double &limitErr, const double *hint)
 
void setFitResultTrees (const RooArgSet *, double *)
 
- Protected Member Functions inherited from FitterAlgoBase
RooFitResult * doFit (RooAbsPdf &pdf, RooAbsData &data, RooRealVar &r, const RooCmdArg &constrain, bool doHesse=true, int ndim=1, bool reuseNLL=false)
 
RooFitResult * doFit (RooAbsPdf &pdf, RooAbsData &data, const RooArgList &rs, const RooCmdArg &constrain, bool doHesse=true, int ndim=1, bool reuseNLL=false)
 
double findCrossing (CascadeMinimizer &minim, RooAbsReal &nll, RooRealVar &r, double level, double rStart, double rBound)
 

Protected Attributes

int currentToy_
 
std::auto_ptr< TFile > fitOut
 
int fitStatus_
 
double * globalObservables_
 
double mu_
 
double nll_bonly_
 
double nll_nll0_
 
double nll_sb_
 
int nToys
 
double * nuisanceParameters_
 
int numbadnll_
 
TTree * t_fit_b_
 
TTree * t_fit_sb_
 
- Protected Attributes inherited from FitterAlgoBase
std::auto_ptr< RooAbsReal > nll
 
- Protected Attributes inherited from LimitAlgo
boost::program_options::options_description options_
 

Static Protected Attributes

static std::string backgroundPdfNames_ = "shapeBkg*"
 
static bool justFit_ = false
 
static bool makePlots_ = false
 
static std::string minos_ = "poi"
 
static std::string name_ = ""
 
static bool noErrors_ = false
 
static std::string out_ = "."
 
static float rebinFactor_ = 1.0
 
static bool reuseParams_ = false
 
static bool saveNormalizations_ = false
 
static std::string signalPdfNames_ = "shapeSig*"
 
- Static Protected Attributes inherited from FitterAlgoBase
static bool do95_ = false
 
static bool keepFailures_ = false
 
static int maxFailedSteps_ = 5
 
static std::string minimizerAlgo_ = "Minuit2"
 
static std::string minimizerAlgoForMinos_ = "Minuit2,simplex"
 
static int minimizerStrategy_ = 1
 
static int minimizerStrategyForMinos_ = 0
 
static float minimizerTolerance_ = 1e-2
 
static float minimizerToleranceForMinos_ = 1e-4
 
static float nllValue_ = std::numeric_limits<float>::quiet_NaN()
 
static float preFitValue_ = 1.0
 
static bool robustFit_ = false
 
static bool saveNLL_ = false
 
static float stepSize_ = 0.1
 

Detailed Description

Do a ML fit of the data with background and signal+background hypothesis and print out diagnostics plots

Author
Giovanni Petrucciani (UCSD)

Definition at line 13 of file MaxLikelihoodFit.h.

Constructor & Destructor Documentation

MaxLikelihoodFit::MaxLikelihoodFit ( )

Definition at line 43 of file MaxLikelihoodFit.cc.

References backgroundPdfNames_, fitStatus_, minos_, mu_, nll_bonly_, nll_nll0_, nll_sb_, nToys, numbadnll_, LimitAlgo::options_, out_, rebinFactor_, and signalPdfNames_.

43  :
44  FitterAlgoBase("MaxLikelihoodFit specific options")
45 {
46  options_.add_options()
47  ("minos", boost::program_options::value<std::string>(&minos_)->default_value(minos_), "Compute MINOS errors for: 'none', 'poi', 'all'")
48  ("out", boost::program_options::value<std::string>(&out_)->default_value(out_), "Directory to put output in")
49  ("plots", "Make plots")
50  ("rebinFactor", boost::program_options::value<float>(&rebinFactor_)->default_value(rebinFactor_), "Rebin by this factor before plotting (does not affect fitting!)")
51  ("signalPdfNames", boost::program_options::value<std::string>(&signalPdfNames_)->default_value(signalPdfNames_), "Names of signal pdfs in plots (separated by ,)")
52  ("backgroundPdfNames", boost::program_options::value<std::string>(&backgroundPdfNames_)->default_value(backgroundPdfNames_), "Names of background pdfs in plots (separated by ',')")
53  ("saveNormalizations", "Save post-fit normalizations of all components of the pdfs")
54  ("justFit", "Just do the S+B fit, don't do the B-only one, don't save output file")
55  ("noErrors", "Don't compute uncertainties on the best fit value")
56  ("initFromBonly", "Use the vlaues of the nuisance parameters from the background only fit as the starting point for the s+b fit")
57  ;
58 
59  // setup a few defaults
61 }
static std::string backgroundPdfNames_
static float rebinFactor_
FitterAlgoBase(const char *title="<FillMe> specific options")
static std::string minos_
static std::string signalPdfNames_
static std::string out_
boost::program_options::options_description options_
Definition: LimitAlgo.h:32
MaxLikelihoodFit::~MaxLikelihoodFit ( )

Definition at line 63 of file MaxLikelihoodFit.cc.

References globalObservables_, and nuisanceParameters_.

63  {
64  // delete the Arrays used to fill the trees;
65  delete globalObservables_;
66  delete nuisanceParameters_;
67 }
double * nuisanceParameters_
double * globalObservables_

Member Function Documentation

void MaxLikelihoodFit::applyOptions ( const boost::program_options::variables_map &  vm)
virtual

Reimplemented from LimitAlgo.

Definition at line 75 of file MaxLikelihoodFit.cc.

References FitterAlgoBase::applyOptionsBase(), justFit_, makePlots_, name_, noErrors_, out_, reuseParams_, and saveNormalizations_.

76 {
77  applyOptionsBase(vm);
78  makePlots_ = vm.count("plots");
79  name_ = vm["name"].defaulted() ? std::string() : vm["name"].as<std::string>();
80  saveNormalizations_ = vm.count("saveNormalizations");
81  justFit_ = vm.count("justFit");
82  noErrors_ = vm.count("noErrors");
83  reuseParams_ = vm.count("initFromBonly");
84  if (justFit_) { out_ = "none"; makePlots_ = false; saveNormalizations_ = false; reuseParams_ = false;}
85  // For now default this to true;
86 }
void applyOptionsBase(const boost::program_options::variables_map &vm)
static bool saveNormalizations_
static bool reuseParams_
static bool justFit_
static bool noErrors_
static bool makePlots_
static std::string out_
static std::string name_
void MaxLikelihoodFit::createFitResultTrees ( const RooStats::ModelConfig &  mc)
protected

Definition at line 402 of file MaxLikelihoodFit.cc.

References a, prof2calltree::count, gather_cfg::cout, fitStatus_, globalObservables_, mu_, name(), nll_bonly_, nll_nll0_, nll_sb_, nuisanceParameters_, numbadnll_, t_fit_b_, and t_fit_sb_.

Referenced by runSpecific().

402  {
403 
404  // Initiate the arrays to store parameters
405 
406  // create TTrees to store fit results:
407  t_fit_b_ = new TTree("tree_fit_b","tree_fit_b");
408  t_fit_sb_ = new TTree("tree_fit_sb","tree_fit_sb");
409 
410  t_fit_b_->Branch("fit_status",&fitStatus_,"fit_status/Int_t");
411  t_fit_sb_->Branch("fit_status",&fitStatus_,"fit_status/Int_t");
412 
413  t_fit_b_->Branch("mu",&mu_,"mu/Double_t");
414  t_fit_sb_->Branch("mu",&mu_,"mu/Double_t");
415 
416  t_fit_b_->Branch("numbadnll",&numbadnll_,"numbadnll/Int_t");
417  t_fit_sb_->Branch("numbadnll",&numbadnll_,"numbadnll/Int_t");
418 
419  t_fit_b_->Branch("nll_min",&nll_bonly_,"nll_min/Double_t");
420  t_fit_sb_->Branch("nll_min",&nll_sb_,"nll_min/Double_t");
421 
422  t_fit_sb_->Branch("nll_nll0",&nll_nll0_,"nll_nll0/Double_t");
423 
424  // fill the maps for the nuisances, and global observables
425  const RooArgSet *cons = mc.GetGlobalObservables();
426  const RooArgSet *nuis = mc.GetNuisanceParameters();
427 
428  globalObservables_ = new double[cons->getSize()];
429  nuisanceParameters_= new double[nuis->getSize()];
430 
431  int count=0;
432  TIterator* iter_c(cons->createIterator());
433  for (TObject *a = iter_c->Next(); a != 0; a = iter_c->Next()) {
434  RooRealVar *rrv = dynamic_cast<RooRealVar *>(a);
435  std::string name = rrv->GetName();
437  t_fit_sb_->Branch(name.c_str(),&(globalObservables_[count]),Form("%s/Double_t",name.c_str()));
438  t_fit_b_->Branch(name.c_str(),&(globalObservables_[count]),Form("%s/Double_t",name.c_str()));
439  count++;
440  }
441 
442  count = 0;
443  TIterator* iter_n(nuis->createIterator());
444  for (TObject *a = iter_n->Next(); a != 0; a = iter_n->Next()) {
445  RooRealVar *rrv = dynamic_cast<RooRealVar *>(a);
446  std::string name = rrv->GetName();
448  t_fit_sb_->Branch(name.c_str(),&(nuisanceParameters_[count])),Form("%s/Double_t",name.c_str());
449  t_fit_b_->Branch(name.c_str(),&(nuisanceParameters_[count]),Form("%s/Double_t",name.c_str()));
450  count++;
451  }
452  std::cout << "Created Branches" <<std::endl;
453  return;
454 }
double * nuisanceParameters_
double a
Definition: hdecay.h:121
tuple cout
Definition: gather_cfg.py:121
double * globalObservables_
virtual const std::string & name() const
void MaxLikelihoodFit::getNormalizations ( RooAbsPdf *  pdf,
const RooArgSet &  obs,
RooArgSet &  out 
)
protected

Definition at line 355 of file MaxLikelihoodFit.cc.

References Clusterizer1DCommons::add(), i, list(), n, and parseEventContent::prod.

Referenced by runSpecific().

355  {
356  RooSimultaneous *sim = dynamic_cast<RooSimultaneous *>(pdf);
357  if (sim != 0) {
358  RooAbsCategoryLValue &cat = const_cast<RooAbsCategoryLValue &>(sim->indexCat());
359  for (int i = 0, n = cat.numBins((const char *)0); i < n; ++i) {
360  cat.setBin(i);
361  RooAbsPdf *pdfi = sim->getPdf(cat.getLabel());
362  if (pdfi) getNormalizations(pdfi, obs, out);
363  }
364  return;
365  }
366  RooProdPdf *prod = dynamic_cast<RooProdPdf *>(pdf);
367  if (prod != 0) {
368  RooArgList list(prod->pdfList());
369  for (int i = 0, n = list.getSize(); i < n; ++i) {
370  RooAbsPdf *pdfi = (RooAbsPdf *) list.at(i);
371  if (pdfi->dependsOn(obs)) getNormalizations(pdfi, obs, out);
372  }
373  return;
374  }
375  RooAddPdf *add = dynamic_cast<RooAddPdf *>(pdf);
376  if (add != 0) {
377  RooArgList list(add->coefList());
378  for (int i = 0, n = list.getSize(); i < n; ++i) {
379  RooAbsReal *coeff = (RooAbsReal *) list.at(i);
380  out.addOwned(*(new RooConstVar(coeff->GetName(), "", coeff->getVal())));
381  }
382  return;
383  }
384 }
int i
Definition: DBlmapReader.cc:9
void add(const std::vector< const T * > &source, std::vector< const T * > &dest)
Definition: sim.h:19
tuple out
Definition: dbtoconf.py:99
void getNormalizations(RooAbsPdf *pdf, const RooArgSet &obs, RooArgSet &out)
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger list("!*","!HLTx*"if it matches 2 triggers or more) will accept the event if all the matching triggers are FAIL.It will reject the event if any of the triggers are PASS or EXCEPTION(this matches the behavior of"!*"before the partial wildcard feature was incorporated).Triggers which are in the READY state are completely ignored.(READY should never be returned since the trigger paths have been run
virtual const std::string& MaxLikelihoodFit::name ( void  ) const
inlinevirtual

Implements LimitAlgo.

Definition at line 16 of file MaxLikelihoodFit.h.

Referenced by createFitResultTrees(), and setFitResultTrees().

16  {
17  static const std::string name("MaxLikelihoodFit");
18  return name;
19  }
virtual const std::string & name() const
bool MaxLikelihoodFit::runSpecific ( RooWorkspace *  w,
RooStats::ModelConfig *  mc_s,
RooStats::ModelConfig *  mc_b,
RooAbsData &  data,
double &  limit,
double &  limitErr,
const double *  hint 
)
protectedvirtual

Implements FitterAlgoBase.

Definition at line 88 of file MaxLikelihoodFit.cc.

References backgroundPdfNames_, alignmentValidation::c1, Combine::commitPoint(), corr, gather_cfg::cout, createFitResultTrees(), currentToy_, data, FitterAlgoBase::do95_, FitterAlgoBase::doFit(), fitOut, fitStatus_, getNormalizations(), globalObservables_, justFit_, utils::makeNuisancePdf(), utils::makePlots(), makePlots_, FitterAlgoBase::minimizerStrategy_, minos_, mu_, name_, FitterAlgoBase::nll, nll_bonly_, nll_nll0_, nll_sb_, noErrors_, nToys, nuisanceParameters_, numbadnll_, out_, RecoTauValidation_cfi::plots, FitterAlgoBase::preFitValue_, alignCSCRings::r, rebinFactor_, reuseParams_, saveNormalizations_, setFitResultTrees(), signalPdfNames_, t_fit_b_, t_fit_sb_, utils::tdrStyle(), and withSystematics.

88  {
89 
90  if (reuseParams_ && minos_!="none"){
91  std::cout << "Cannot reuse b-only fit params when running minos. Parameters will be reset when running S+B fit"<<std::endl;
92  reuseParams_=false;
93  }
94 
95  if (!justFit_ && out_ != "none"){
96  if (currentToy_ < 1){
97  fitOut.reset(TFile::Open((out_+"/mlfit"+name_+".root").c_str(), "RECREATE"));
98  createFitResultTrees(*mc_s);
99  }
100  }
101 
102  RooRealVar *r = dynamic_cast<RooRealVar *>(mc_s->GetParametersOfInterest()->first());
103 
104  TCanvas *c1 = 0;
105  if (makePlots_) {
106  utils::tdrStyle();
107  c1 = new TCanvas("c1","c1");
108  }
109 
110  // Make pre-plots before the fit
111  r->setVal(preFitValue_);
112  if (makePlots_) {
113  std::vector<RooPlot *> plots = utils::makePlots(*mc_s->GetPdf(), data, signalPdfNames_.c_str(), backgroundPdfNames_.c_str(), rebinFactor_);
114  for (std::vector<RooPlot *>::iterator it = plots.begin(), ed = plots.end(); it != ed; ++it) {
115  (*it)->Draw();
116  c1->Print((out_+"/"+(*it)->GetName()+"_prefit.png").c_str());
117  if (fitOut.get() && currentToy_< 1) fitOut->WriteTObject(*it, (std::string((*it)->GetName())+"_prefit").c_str());
118  }
119  }
120 
121 
122  // Determine pre-fit values of nuisance parameters
123  if (currentToy_ < 1){
124  const RooArgSet *nuis = mc_s->GetNuisanceParameters();
125  const RooArgSet *globalObs = mc_s->GetGlobalObservables();
126  if (!justFit_ && nuis && globalObs ) {
127  std::auto_ptr<RooAbsPdf> nuisancePdf(utils::makeNuisancePdf(*mc_s));
128  std::auto_ptr<RooDataSet> globalData(new RooDataSet("globalData","globalData", *globalObs));
129  globalData->add(*globalObs);
130  RooFitResult *res_prefit = 0;
131  {
132  CloseCoutSentry sentry(verbose < 2);
133  res_prefit = nuisancePdf->fitTo(*globalData,
134  RooFit::Save(1),
135  RooFit::Minimizer(ROOT::Math::MinimizerOptions::DefaultMinimizerType().c_str(), ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo().c_str()),
136  RooFit::Strategy(minimizerStrategy_),
137  RooFit::Minos(minos_ == "all")
138  );
139  }
140  if (fitOut.get() ) fitOut->WriteTObject(res_prefit, "nuisances_prefit_res");
141  if (fitOut.get() ) fitOut->WriteTObject(nuis->snapshot(), "nuisances_prefit");
142 
143  nuisancePdf.reset();
144  globalData.reset();
145  delete res_prefit;
146 
147  } else if (nuis) {
148  if (fitOut.get() ) fitOut->WriteTObject(nuis->snapshot(), "nuisances_prefit");
149  }
150  }
151 
152  RooFitResult *res_b = 0, *res_s = 0;
153  const RooCmdArg &constCmdArg_s = withSystematics ? RooFit::Constrain(*mc_s->GetNuisanceParameters()) : RooFit::NumCPU(1); // use something dummy
154  const RooCmdArg &minosCmdArg = minos_ == "poi" ? RooFit::Minos(*mc_s->GetParametersOfInterest()) : RooFit::Minos(minos_ != "none");
155  w->loadSnapshot("clean");
156  r->setVal(0.0); r->setConstant(true);
157 
158  // Setup Nll before calling fits;
159  if (currentToy_<1) nll.reset(mc_s->GetPdf()->createNLL(data,constCmdArg_s,RooFit::Extended(mc_s->GetPdf()->canBeExtended())));
160  // Get the nll value on the prefit
161  double nll0 = nll->getVal();
162 
163  if (justFit_) {
164  // skip b-only fit
165  } else if (minos_ != "all") {
166  RooArgList minos;
167  res_b = doFit(*mc_s->GetPdf(), data, minos, constCmdArg_s, /*hesse=*/true,/*reuseNLL*/ true);
168  nll_bonly_=nll->getVal()-nll0;
169  } else {
170  CloseCoutSentry sentry(verbose < 2);
171  res_b = mc_s->GetPdf()->fitTo(data,
172  RooFit::Save(1),
173  RooFit::Minimizer(ROOT::Math::MinimizerOptions::DefaultMinimizerType().c_str(), ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo().c_str()),
174  RooFit::Strategy(minimizerStrategy_),
175  RooFit::Extended(mc_s->GetPdf()->canBeExtended()),
176  constCmdArg_s, minosCmdArg
177  );
178  if (res_b) nll_bonly_ = nll->getVal() - nll0;
179 
180  }
181 
182  if (res_b) {
183  if (verbose > 1) res_b->Print("V");
184  if (fitOut.get()) {
185  if (currentToy_< 1) fitOut->WriteTObject(res_b,"fit_b");
186  setFitResultTrees(mc_s->GetNuisanceParameters(),nuisanceParameters_);
187  setFitResultTrees(mc_s->GetGlobalObservables(),globalObservables_);
188  fitStatus_ = res_b->status();
189  }
190  numbadnll_=res_b->numInvalidNLL();
191 
192  if (makePlots_) {
193  std::vector<RooPlot *> plots = utils::makePlots(*mc_b->GetPdf(), data, signalPdfNames_.c_str(), backgroundPdfNames_.c_str(), rebinFactor_);
194  for (std::vector<RooPlot *>::iterator it = plots.begin(), ed = plots.end(); it != ed; ++it) {
195  c1->cd(); (*it)->Draw();
196  c1->Print((out_+"/"+(*it)->GetName()+"_fit_b.png").c_str());
197  if (fitOut.get() && currentToy_< 1) fitOut->WriteTObject(*it, (std::string((*it)->GetName())+"_fit_b").c_str());
198  }
199  }
200 
201  if (saveNormalizations_ && currentToy_<1) {
202  RooArgSet *norms = new RooArgSet();
203  norms->setName("norm_fit_b");
204  getNormalizations(mc_s->GetPdf(), *mc_s->GetObservables(), *norms);
205  if (fitOut.get()) fitOut->WriteTObject(norms, "norm_fit_b");
206  delete norms;
207  }
208 
209  if (makePlots_ && currentToy_<1) {
210  TH2 *corr = res_b->correlationHist();
211  c1->SetLeftMargin(0.25); c1->SetBottomMargin(0.25);
212  corr->SetTitle("Correlation matrix of fit parameters");
213  gStyle->SetPaintTextFormat(res_b->floatParsFinal().getSize() > 10 ? ".1f" : ".2f");
214  gStyle->SetOptStat(0);
215  corr->SetMarkerSize(res_b->floatParsFinal().getSize() > 10 ? 2 : 1);
216  corr->Draw("COLZ TEXT");
217  c1->Print((out_+"/covariance_fit_b.png").c_str());
218  c1->SetLeftMargin(0.16); c1->SetBottomMargin(0.13);
219  if (fitOut.get()) fitOut->WriteTObject(corr, "covariance_fit_b");
220  }
221  }
222  else {
223  fitStatus_=-1;
224  numbadnll_=-1;
225  }
226  mu_=r->getVal();
227  if (t_fit_b_) t_fit_b_->Fill();
228  // no longer need res_b
229  delete res_b;
230 
231  if (!reuseParams_) w->loadSnapshot("clean"); // Reset, also ensures nll_prefit is same in call to doFit for b and s+b
232  r->setVal(preFitValue_); r->setConstant(false);
233  if (minos_ != "all") {
234  RooArgList minos; if (minos_ == "poi") minos.add(*r);
235  res_s = doFit(*mc_s->GetPdf(), data, minos, constCmdArg_s, /*hesse=*/!noErrors_,/*reuseNLL*/ true);
236  nll_sb_ = nll->getVal()-nll0;
237  } else {
238  CloseCoutSentry sentry(verbose < 2);
239  res_s = mc_s->GetPdf()->fitTo(data,
240  RooFit::Save(1),
241  RooFit::Minimizer(ROOT::Math::MinimizerOptions::DefaultMinimizerType().c_str(), ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo().c_str()),
242  RooFit::Strategy(minimizerStrategy_),
243  RooFit::Extended(mc_s->GetPdf()->canBeExtended()),
244  constCmdArg_s, minosCmdArg
245  );
246  if (res_s) nll_sb_= nll->getVal()-nll0;
247 
248  }
249  if (res_s) {
250  limit = r->getVal();
251  limitErr = r->getError();
252  if (verbose > 1) res_s->Print("V");
253  if (fitOut.get()){
254  if (currentToy_<1) fitOut->WriteTObject(res_s, "fit_s");
255 
256  setFitResultTrees(mc_s->GetNuisanceParameters(),nuisanceParameters_);
257  setFitResultTrees(mc_s->GetGlobalObservables(),globalObservables_);
258  fitStatus_ = res_s->status();
259  numbadnll_ = res_s->numInvalidNLL();
260 
261  // Additionally store the nll_sb - nll_bonly (=0.5*q0)
262  nll_nll0_ = nll_sb_ - nll_bonly_;
263  }
264 
265  if (makePlots_) {
266  std::vector<RooPlot *> plots = utils::makePlots(*mc_s->GetPdf(), data, signalPdfNames_.c_str(), backgroundPdfNames_.c_str(), rebinFactor_);
267  for (std::vector<RooPlot *>::iterator it = plots.begin(), ed = plots.end(); it != ed; ++it) {
268  c1->cd(); (*it)->Draw();
269  c1->Print((out_+"/"+(*it)->GetName()+"_fit_s.png").c_str());
270  if (fitOut.get() && currentToy_< 1) fitOut->WriteTObject(*it, (std::string((*it)->GetName())+"_fit_s").c_str());
271  }
272  }
273 
274  if (saveNormalizations_&& currentToy_< 1) {
275  RooArgSet *norms = new RooArgSet();
276  norms->setName("norm_fit_s");
277  getNormalizations(mc_s->GetPdf(), *mc_s->GetObservables(), *norms);
278  if (fitOut.get() ) fitOut->WriteTObject(norms, "norm_fit_s");
279  delete norms;
280  }
281 
282  if (makePlots_&& currentToy_< 1) {
283  TH2 *corr = res_s->correlationHist();
284  c1->SetLeftMargin(0.25); c1->SetBottomMargin(0.25);
285  corr->SetTitle("Correlation matrix of fit parameters");
286  gStyle->SetPaintTextFormat(res_s->floatParsFinal().getSize() > 10 ? ".1f" : ".2f");
287  gStyle->SetOptStat(0);
288  corr->SetMarkerSize(res_s->floatParsFinal().getSize() > 10 ? 2 : 1);
289  corr->Draw("COLZ TEXT");
290  c1->Print((out_+"/covariance_fit_s.png").c_str());
291  c1->SetLeftMargin(0.16); c1->SetBottomMargin(0.13);
292  if (fitOut.get() ) fitOut->WriteTObject(corr, "covariance_fit_s");
293  }
294  } else {
295  fitStatus_=-1;
296  numbadnll_=-1;
297  nll_nll0_ = -1;
298  }
299  mu_=r->getVal();
300  if (t_fit_sb_) t_fit_sb_->Fill();
301 
302  if (res_s) {
303  RooRealVar *rf = dynamic_cast<RooRealVar*>(res_s->floatParsFinal().find(r->GetName()));
304  double bestFitVal = rf->getVal();
305 
306  double hiErr = +(rf->hasRange("err68") ? rf->getMax("err68") - bestFitVal : rf->getAsymErrorHi());
307  double loErr = -(rf->hasRange("err68") ? rf->getMin("err68") - bestFitVal : rf->getAsymErrorLo());
308  double maxError = std::max<double>(std::max<double>(hiErr, loErr), rf->getError());
309 
310  if (fabs(hiErr) < 0.001*maxError) hiErr = -bestFitVal + rf->getMax();
311  if (fabs(loErr) < 0.001*maxError) loErr = +bestFitVal - rf->getMin();
312 
313  double hiErr95 = +(do95_ && rf->hasRange("err95") ? rf->getMax("err95") - bestFitVal : 0);
314  double loErr95 = -(do95_ && rf->hasRange("err95") ? rf->getMin("err95") - bestFitVal : 0);
315 
316  limit = bestFitVal; limitErr = 0;
317  if (!noErrors_) Combine::commitPoint(/*expected=*/true, /*quantile=*/0.5);
318  limit = bestFitVal - loErr; limitErr = 0;
319  if (!noErrors_) Combine::commitPoint(/*expected=*/true, /*quantile=*/0.16);
320  limit = bestFitVal + hiErr; limitErr = 0;
321  if (!noErrors_) Combine::commitPoint(/*expected=*/true, /*quantile=*/0.84);
322  if (do95_ && rf->hasRange("err95") && !noErrors_) {
323  limit = rf->getMax("err95"); Combine::commitPoint(/*expected=*/true, /*quantile=*/0.975);
324  limit = rf->getMin("err95"); Combine::commitPoint(/*expected=*/true, /*quantile=*/0.025);
325  }
326 
327  limit = bestFitVal;
328  limitErr = maxError;
329  std::cout << "\n --- MaxLikelihoodFit ---" << std::endl;
330  std::cout << "Best fit " << r->GetName() << ": " << rf->getVal() << " "<< -loErr << "/+" << +hiErr << " (68% CL)" << std::endl;
331  if (do95_) {
332  std::cout << " " << r->GetName() << ": " << rf->getVal() << " "<< -loErr95 << "/+" << +hiErr95 << " (95% CL)" << std::endl;
333  }
334  } else {
335  std::cout << "\n --- MaxLikelihoodFit ---" << std::endl;
336  std::cout << "Fit failed." << std::endl;
337  }
338 
339  if (currentToy_==nToys-1 || nToys==0 ) {
340 
341  if (fitOut.get()) {
342  fitOut->cd();
343  t_fit_sb_->Write(); t_fit_b_->Write();
344  fitOut.release()->Close();
345  }
346 
347  }
348  bool fitreturn = (res_s!=0);
349  delete res_s;
350 
351  std::cout << "nll S+B -> "<<nll_sb_ << " nll B -> " << nll_bonly_ <<std::endl;
352  return fitreturn;
353 }
void createFitResultTrees(const RooStats::ModelConfig &)
double * nuisanceParameters_
void setFitResultTrees(const RooArgSet *, double *)
std::vector< RooPlot * > makePlots(const RooAbsPdf &pdf, const RooAbsData &data, const char *signalSel=0, const char *backgroundSel=0, float rebinFactor=1.0)
make plots, if possible
Definition: utils.cc:403
static bool saveNormalizations_
void tdrStyle()
set style for plots
Definition: tdrstyle.cc:4
bool withSystematics
Definition: Combine.cc:68
std::auto_ptr< TFile > fitOut
RooAbsPdf * makeNuisancePdf(RooStats::ModelConfig &model, const char *name="nuisancePdf")
Definition: utils.cc:200
static bool do95_
static bool reuseParams_
static bool justFit_
static std::string backgroundPdfNames_
RooFitResult * doFit(RooAbsPdf &pdf, RooAbsData &data, RooRealVar &r, const RooCmdArg &constrain, bool doHesse=true, int ndim=1, bool reuseNLL=false)
static float rebinFactor_
static void commitPoint(bool expected, float quantile)
Save a point into the output tree. Usually if expected = false, quantile should be set to -1 (except ...
Definition: Combine.cc:557
static float preFitValue_
std::auto_ptr< RooAbsReal > nll
static bool noErrors_
static int minimizerStrategy_
JetCorrectorParameters corr
Definition: classes.h:9
static std::string minos_
static bool makePlots_
static std::string signalPdfNames_
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
static std::string out_
static std::string name_
tuple cout
Definition: gather_cfg.py:121
double * globalObservables_
void getNormalizations(RooAbsPdf *pdf, const RooArgSet &obs, RooArgSet &out)
T w() const
void MaxLikelihoodFit::setFitResultTrees ( const RooArgSet *  args,
double *  vals 
)
protected

Definition at line 387 of file MaxLikelihoodFit.cc.

References a, prof2calltree::count, and name().

Referenced by runSpecific().

387  {
388 
389  TIterator* iter(args->createIterator());
390  int count=0;
391 
392  for (TObject *a = iter->Next(); a != 0; a = iter->Next()) {
393  RooRealVar *rrv = dynamic_cast<RooRealVar *>(a);
394  std::string name = rrv->GetName();
395  vals[count]=rrv->getVal();
396  count++;
397  }
398  delete iter;
399  return;
400 }
dictionary args
double a
Definition: hdecay.h:121
virtual const std::string & name() const
void MaxLikelihoodFit::setNToys ( const int  iToy)
virtual

Reimplemented from LimitAlgo.

Definition at line 72 of file MaxLikelihoodFit.cc.

References nToys.

72  {
73  nToys = iToy;
74 }
void MaxLikelihoodFit::setToyNumber ( const int  iToy)
virtual

Reimplemented from LimitAlgo.

Definition at line 69 of file MaxLikelihoodFit.cc.

References currentToy_.

69  {
70  currentToy_ = iToy;
71 }

Member Data Documentation

std::string MaxLikelihoodFit::backgroundPdfNames_ = "shapeBkg*"
staticprotected

Definition at line 36 of file MaxLikelihoodFit.h.

Referenced by MaxLikelihoodFit(), and runSpecific().

int MaxLikelihoodFit::currentToy_
protected

Definition at line 39 of file MaxLikelihoodFit.h.

Referenced by runSpecific(), and setToyNumber().

std::auto_ptr<TFile> MaxLikelihoodFit::fitOut
protected

Definition at line 42 of file MaxLikelihoodFit.h.

Referenced by runSpecific().

int MaxLikelihoodFit::fitStatus_
protected

Definition at line 40 of file MaxLikelihoodFit.h.

Referenced by createFitResultTrees(), MaxLikelihoodFit(), and runSpecific().

double* MaxLikelihoodFit::globalObservables_
protected

Definition at line 43 of file MaxLikelihoodFit.h.

Referenced by createFitResultTrees(), runSpecific(), and ~MaxLikelihoodFit().

bool MaxLikelihoodFit::justFit_ = false
staticprotected

Definition at line 32 of file MaxLikelihoodFit.h.

Referenced by applyOptions(), and runSpecific().

bool MaxLikelihoodFit::makePlots_ = false
staticprotected

Definition at line 34 of file MaxLikelihoodFit.h.

Referenced by applyOptions(), and runSpecific().

std::string MaxLikelihoodFit::minos_ = "poi"
staticprotected

Definition at line 30 of file MaxLikelihoodFit.h.

Referenced by MaxLikelihoodFit(), and runSpecific().

double MaxLikelihoodFit::mu_
protected

Definition at line 41 of file MaxLikelihoodFit.h.

Referenced by createFitResultTrees(), MaxLikelihoodFit(), and runSpecific().

std::string MaxLikelihoodFit::name_ = ""
staticprotected

Definition at line 28 of file MaxLikelihoodFit.h.

Referenced by applyOptions(), and runSpecific().

double MaxLikelihoodFit::nll_bonly_
protected

Definition at line 41 of file MaxLikelihoodFit.h.

Referenced by createFitResultTrees(), MaxLikelihoodFit(), and runSpecific().

double MaxLikelihoodFit::nll_nll0_
protected

Definition at line 41 of file MaxLikelihoodFit.h.

Referenced by createFitResultTrees(), MaxLikelihoodFit(), and runSpecific().

double MaxLikelihoodFit::nll_sb_
protected

Definition at line 41 of file MaxLikelihoodFit.h.

Referenced by createFitResultTrees(), MaxLikelihoodFit(), and runSpecific().

bool MaxLikelihoodFit::noErrors_ = false
staticprotected

Definition at line 32 of file MaxLikelihoodFit.h.

Referenced by applyOptions(), and runSpecific().

int MaxLikelihoodFit::nToys
protected

Definition at line 39 of file MaxLikelihoodFit.h.

Referenced by MaxLikelihoodFit(), runSpecific(), and setNToys().

double* MaxLikelihoodFit::nuisanceParameters_
protected

Definition at line 44 of file MaxLikelihoodFit.h.

Referenced by createFitResultTrees(), runSpecific(), and ~MaxLikelihoodFit().

int MaxLikelihoodFit::numbadnll_
protected

Definition at line 40 of file MaxLikelihoodFit.h.

Referenced by createFitResultTrees(), MaxLikelihoodFit(), and runSpecific().

std::string MaxLikelihoodFit::out_ = "."
staticprotected

Definition at line 33 of file MaxLikelihoodFit.h.

Referenced by applyOptions(), MaxLikelihoodFit(), and runSpecific().

float MaxLikelihoodFit::rebinFactor_ = 1.0
staticprotected

Definition at line 35 of file MaxLikelihoodFit.h.

Referenced by MaxLikelihoodFit(), and runSpecific().

bool MaxLikelihoodFit::reuseParams_ = false
staticprotected

Definition at line 38 of file MaxLikelihoodFit.h.

Referenced by applyOptions(), and runSpecific().

bool MaxLikelihoodFit::saveNormalizations_ = false
staticprotected

Definition at line 37 of file MaxLikelihoodFit.h.

Referenced by applyOptions(), and runSpecific().

std::string MaxLikelihoodFit::signalPdfNames_ = "shapeSig*"
staticprotected

Definition at line 36 of file MaxLikelihoodFit.h.

Referenced by MaxLikelihoodFit(), and runSpecific().

TTree* MaxLikelihoodFit::t_fit_b_
protected

Definition at line 46 of file MaxLikelihoodFit.h.

Referenced by createFitResultTrees(), and runSpecific().

TTree * MaxLikelihoodFit::t_fit_sb_
protected

Definition at line 46 of file MaxLikelihoodFit.h.

Referenced by createFitResultTrees(), and runSpecific().