Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
RooMCStudy.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitCore *
4 * @(#)root/roofitcore:$Id$
5 * Authors: *
6 * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
7 * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
8 * *
9 * Copyright (c) 2000-2005, Regents of the University of California *
10 * and Stanford University. All rights reserved. *
11 * *
12 * Redistribution and use in source and binary forms, *
13 * with or without modification, are permitted according to the terms *
14 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
15 *****************************************************************************/
16
17/**
18\file RooMCStudy.cxx
19\class RooMCStudy
20\ingroup Roofitcore
21
22Helper class to facilitate Monte Carlo studies
23such as 'goodness-of-fit' studies, that involve fitting a PDF
24to multiple toy Monte Carlo sets. These may be generated from either same PDF
25or from a different PDF with similar parameters.
26
27Given a fit and a generator PDF (they might be identical), RooMCStudy can produce
28toyMC samples and/or fit these.
29It accumulates the post-fit parameters of each iteration in a dataset. These can be
30retrieved using fitParams() or fitParDataSet(). This dataset additionally contains the
31variables
32- NLL: The value of the negative log-likelihood for each run.
33- ngen: The number of events generated for each run.
34
35Additional plotting routines simplify the task of plotting
36the distribution of the minimized likelihood, the fitted parameter values,
37fitted error and pull distribution.
38
39RooMCStudy provides the option to insert add-in modules
40that modify the generate-and-fit cycle and allow to perform
41extra steps in the cycle. Output of these modules can be stored
42alongside the fit results in the aggregate results dataset.
43These study modules should derive from the class RooAbsMCStudyModule.
44
45Check the RooFit tutorials
46- rf801_mcstudy.C
47- rf802_mcstudy_addons.C
48- rf803_mcstudy_addons2.C
49- rf804_mcstudy_constr.C
50for usage examples.
51**/
52
53
54#include <RooMCStudy.h>
55
56#include <RooAbsMCStudyModule.h>
57#include <RooAbsPdf.h>
58#include <RooArgList.h>
59#include <RooCmdConfig.h>
60#include <RooDataHist.h>
61#include <RooDataSet.h>
62#include <RooErrorVar.h>
63#include <RooFitResult.h>
64#include <RooFormulaVar.h>
65#include <RooGenContext.h>
66#include <RooGlobalFunc.h>
67#include <RooMsgService.h>
68#include <RooPlot.h>
69#include <RooProdPdf.h>
70#include <RooPullVar.h>
71#include <RooRandom.h>
72#include <RooRealVar.h>
73#include <RooWorkspace.h>
74
75#include <snprintf.h>
76#include <iostream>
77
79
80/**
81Construct Monte Carlo Study Manager. This class automates generating data from a given PDF,
82fitting the PDF to data and accumulating the fit statistics.
83
84\param[in] model The PDF to be studied
85\param[in] observables The variables of the PDF to be considered observables
86\param[in] arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8 Optional arguments according to table below.
87
88<table>
89<tr><th> Optional arguments <th>
90<tr><td> Silence() <td> Suppress all RooFit messages during running below PROGRESS level
91<tr><td> FitModel(const RooAbsPdf&) <td> The PDF for fitting if it is different from the PDF for generating.
92<tr><td> ConditionalObservables(const RooArgSet& set) <td> The set of observables that the PDF should _not_ be normalized over
93<tr><td> Binned(bool flag) <td> Bin the dataset before fitting it. Speeds up fitting of large data samples
94<tr><td> FitOptions(....) <td> Options to be used for fitting. All named arguments inside FitOptions() are passed to RooAbsPdf::fitTo().
95 `Save()` is especially interesting to be able to retrieve fit results of each run using fitResult().
96<tr><td> Verbose(bool flag) <td> Activate informational messages in event generation phase
97<tr><td> Extended(bool flag) <td> Determine number of events for each sample anew from a Poisson distribution
98<tr><td> Constrain(const RooArgSet& pars) <td> Apply internal constraints on given parameters in fit and sample constrained parameter values from constraint p.d.f for each toy.
99<tr><td> ExternalConstraints(const RooArgSet& ) <td> Apply internal constraints on given parameters in fit and sample constrained parameter values from constraint p.d.f for each toy.
100<tr><td> ProtoData(const RooDataSet&, bool randOrder)
101 <td> Prototype data for the event generation. If the randOrder flag is set, the order of the dataset will be re-randomized for each generation
102 cycle to protect against systematic biases if the number of generated events does not exactly match the number of events in the prototype dataset
103 at the cost of reduced precision with mu equal to the specified number of events
104</table>
105*/
106RooMCStudy::RooMCStudy(const RooAbsPdf& model, const RooArgSet& observables,
107 const RooCmdArg& arg1, const RooCmdArg& arg2,
108 const RooCmdArg& arg3,const RooCmdArg& arg4,const RooCmdArg& arg5,
109 const RooCmdArg& arg6,const RooCmdArg& arg7,const RooCmdArg& arg8) : TNamed("mcstudy","mcstudy")
110
111{
112 // Stuff all arguments in a list
113 RooLinkedList cmdList;
114 cmdList.Add(const_cast<RooCmdArg*>(&arg1)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg2)) ;
115 cmdList.Add(const_cast<RooCmdArg*>(&arg3)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg4)) ;
116 cmdList.Add(const_cast<RooCmdArg*>(&arg5)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg6)) ;
117 cmdList.Add(const_cast<RooCmdArg*>(&arg7)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg8)) ;
118
119 // Select the pdf-specific commands
120 RooCmdConfig pc("RooMCStudy::RooMCStudy(" + std::string(model.GetName()) + ")");
121
122 pc.defineObject("fitModel","FitModel",0,nullptr) ;
123 pc.defineSet("condObs","ProjectedObservables",0,nullptr) ;
124 pc.defineObject("protoData","PrototypeData",0,nullptr) ;
125 pc.defineSet("cPars","Constrain",0,nullptr) ;
126 pc.defineSet("extCons","ExternalConstraints",0,nullptr) ;
127 pc.defineInt("silence","Silence",0,0) ;
128 pc.defineInt("randProtoData","PrototypeData",0,0) ;
129 pc.defineInt("verboseGen","Verbose",0,0) ;
130 pc.defineInt("extendedGen","Extended",0,0) ;
131 pc.defineInt("binGenData","Binned",0,0) ;
132 pc.defineInt("dummy","FitOptArgs",0,0) ;
133
134 // Process and check varargs
135 pc.process(cmdList) ;
136 if (!pc.ok(true)) {
137 // WVE do something here
138 throw std::string("RooMCStudy::RooMCStudy() Error in parsing arguments passed to constructor") ;
139 return ;
140 }
141
142 // Save fit command options
143 if (pc.hasProcessed("FitOptArgs")) {
144 RooCmdArg* fitOptArg = static_cast<RooCmdArg*>(cmdList.FindObject("FitOptArgs")) ;
145 for (int i=0 ; i<fitOptArg->subArgs().GetSize() ;i++) {
146 _fitOptList.Add(new RooCmdArg(static_cast<RooCmdArg&>(*fitOptArg->subArgs().At(i)))) ;
147 }
148 }
149
150 // Decode command line arguments
151 _silence = pc.getInt("silence") ;
152 _verboseGen = pc.getInt("verboseGen") ;
153 _extendedGen = pc.getInt("extendedGen") ;
154 _binGenData = pc.getInt("binGenData") ;
155 _randProto = pc.getInt("randProtoData") ;
156
157 // Process constraints specifications
158 const RooArgSet* cParsTmp = pc.getSet("cPars") ;
159 const RooArgSet* extCons = pc.getSet("extCons") ;
160
161 RooArgSet* cPars = new RooArgSet ;
162 if (cParsTmp) {
163 cPars->add(*cParsTmp) ;
164 }
165
166 // If constraints are specified, add to fit options
167 if (cPars) {
169 }
170 if (extCons) {
172 }
173
174 // Make list of all constraints
175 RooArgSet allConstraints ;
176 RooArgSet consPars ;
177 if (cPars) {
178 if (std::unique_ptr<RooArgSet> constraints{model.getAllConstraints(observables,*cPars,true)}) {
179 allConstraints.add(*constraints) ;
180 }
181 }
182
183 // Construct constraint p.d.f
184 if (!allConstraints.empty()) {
185 _constrPdf = std::make_unique<RooProdPdf>("mcs_constr_prod","RooMCStudy constraints product",allConstraints);
186
187 if (cPars) {
188 consPars.add(*cPars) ;
189 } else {
190 RooArgSet params;
191 model.getParameters(&observables, params);
192 RooArgSet cparams;
193 _constrPdf->getObservables(&params, cparams);
194 consPars.add(cparams) ;
195 }
196 _constrGenContext.reset(_constrPdf->genContext(consPars,nullptr,nullptr,_verboseGen));
197
198 _perExptGenParams = true ;
199
200 coutI(Generation) << "RooMCStudy::RooMCStudy: INFO have pdf with constraints, will generate parameters from constraint pdf for each experiment" << std::endl ;
201 }
202
203
204 // Extract generator and fit models
205 _genModel = const_cast<RooAbsPdf*>(&model) ;
206 RooAbsPdf* fitModel = static_cast<RooAbsPdf*>(pc.getObject("fitModel",nullptr)) ;
207 _fitModel = fitModel ? fitModel : _genModel ;
208
209 // Extract conditional observables and prototype data
210 _genProtoData = static_cast<RooDataSet*>(pc.getObject("protoData",nullptr)) ;
211 if (auto condObs = pc.getSet("condObs",nullptr)) {
212 _projDeps.add(*condObs);
213 }
214
215 _dependents.add(observables) ;
216
218 _canAddFitResults = true ;
219
221 oocoutW(_fitModel,Generation) << "RooMCStudy::RooMCStudy: WARNING Using generator option 'e' (Poisson distribution of #events) together " << std::endl
222 << " with a prototype dataset implies incomplete sampling or oversampling of proto data." << std::endl
223 << " Use option \"r\" to randomize prototype dataset order and thus to randomize" << std::endl
224 << " the set of over/undersampled prototype events for each generation cycle." << std::endl ;
225 }
226
228 if (!_binGenData) {
230 _genContext->attach(_genParams) ;
231 }
232
234
235 // Store list of parameters and save initial values separately
238
240
241 // Place holder for NLL
242 _nllVar = std::make_unique<RooRealVar>("NLL","-log(Likelihood)",0);
243
244 // Place holder for number of generated events
245 _ngenVar = std::make_unique<RooRealVar>("ngen","number of generated events",0);
246
247 // Create data set containing parameter values, errors and pulls
248 RooArgSet tmp2(_fitParams) ;
249 tmp2.add(*_nllVar) ;
250 tmp2.add(*_ngenVar) ;
251
252 // Mark all variable to store their errors in the dataset
253 tmp2.setAttribAll("StoreError",true) ;
254 tmp2.setAttribAll("StoreAsymError",true) ;
255 std::string fpdName;
256 if (_fitModel==_genModel) {
257 fpdName = "fitParData_" + std::string(_fitModel->GetName());
258 } else {
259 fpdName= "fitParData_" + std::string(_fitModel->GetName()) + "_" + std::string(_genModel->GetName());
260 }
261
262 _fitParData = std::make_unique<RooDataSet>(fpdName,"Fit Parameters DataSet",tmp2);
263 tmp2.setAttribAll("StoreError",false) ;
264 tmp2.setAttribAll("StoreAsymError",false) ;
265
266 if (_perExptGenParams) {
267 _genParData = std::make_unique<RooDataSet>("genParData","Generated Parameters dataset",_genParams);
268 }
269
270 // Append proto variables to allDependents
271 if (_genProtoData) {
273 }
274
275 // Call module initializers
276 for (auto iter=_modList.begin() ; iter!= _modList.end() ; ++iter) {
277 bool ok = (*iter)->doInitializeInstance(*this) ;
278 if (!ok) {
279 oocoutE(_fitModel,Generation) << "RooMCStudy::ctor: removing study module " << (*iter)->GetName() << " from analysis chain because initialization failed" << std::endl ;
280 iter = _modList.erase(iter) ;
281 }
282 }
283
284}
285
286
287////////////////////////////////////////////////////////////////////////////////
288
290{
293}
294
295
296
297////////////////////////////////////////////////////////////////////////////////
298/// Insert given RooMCStudy add-on module to the processing chain
299/// of this MCStudy object
300
302{
303 module.doInitializeInstance(*this) ;
304 _modList.push_back(&module) ;
305}
306
307
308
309////////////////////////////////////////////////////////////////////////////////
310/// Run engine method. Generate and/or fit, according to flags, 'nSamples' samples of 'nEvtPerSample' events.
311/// If keepGenData is set, all generated data sets will be kept in memory and can be accessed
312/// later via genData().
313///
314/// When generating, data sets will be written out in ascii form if the pattern string is supplied
315/// The pattern, which is a template for snprintf, should look something like "data/toymc_%04d.dat"
316/// and should contain one integer field that encodes the sample serial number.
317///
318/// When fitting only, data sets may optionally be read from ascii files, using the same file
319/// pattern.
320///
321
322bool RooMCStudy::run(bool doGenerate, bool DoFit, Int_t nSamples, Int_t nEvtPerSample, bool keepGenData, const char* asciiFilePat)
323{
325 if (_silence) {
328 }
329
330 for (RooAbsMCStudyModule *mod : _modList) {
331 mod->initializeRun(nSamples) ;
332 }
333
334 int prescale = nSamples>100 ? int(nSamples/100) : 1 ;
335
336 while(nSamples--) {
337
338 if (nSamples%prescale==0) {
339 oocoutP(_fitModel,Generation) << "RooMCStudy::run: " ;
340 if (doGenerate) ooccoutI(_fitModel,Generation) << "Generating " ;
341 if (doGenerate && DoFit) ooccoutI(_fitModel,Generation) << "and " ;
342 if (DoFit) ooccoutI(_fitModel,Generation) << "fitting " ;
343 ooccoutP(_fitModel,Generation) << "sample " << nSamples << std::endl ;
344 }
345
346 std::unique_ptr<RooAbsData> ownedGenSample;
347 _genSample = nullptr;
348 bool existingData = false ;
349 if (doGenerate) {
350 // Generate sample
351 int nEvt(nEvtPerSample) ;
352
353 // Reset generator parameters to initial values
355
356 // If constraints are present, sample generator values from constraints
357 if (_constrPdf) {
358 _genParams.assign(*std::unique_ptr<RooDataSet>{_constrGenContext->generate(1)}->get());
359 }
360
361 // Save generated parameters if required
362 if (_genParData) {
363 _genParData->add(_genParams) ;
364 }
365
366 // Call module before-generation hook
367 for (RooAbsMCStudyModule *mod : _modList) {
368 mod->processBeforeGen(nSamples) ;
369 }
370
371 if (_binGenData) {
372
373 // Calculate the number of (extended) events for this run
374 if (_extendedGen) {
376 nEvt = RooRandom::randomGenerator()->Poisson(nEvtPerSample==0?_nExpGen:nEvtPerSample) ;
377 }
378
379 // Binned generation
380 ownedGenSample = std::unique_ptr<RooDataHist>{_genModel->generateBinned(_dependents,nEvt)};
381
382 } else {
383
384 // Calculate the number of (extended) events for this run
385 if (_extendedGen) {
387 nEvt = RooRandom::randomGenerator()->Poisson(nEvtPerSample==0?_nExpGen:nEvtPerSample) ;
388 }
389
390 // Optional randomization of protodata for this run
392 oocoutI(_fitModel,Generation) << "RooMCStudy: (Re)randomizing event order in prototype dataset (Nevt=" << nEvt << ")" << std::endl ;
394 _genContext->setProtoDataOrder(newOrder) ;
395 delete[] newOrder ;
396 }
397
398 // Actual generation of events
399 if (nEvt>0) {
400 ownedGenSample = std::unique_ptr<RooAbsData>{_genContext->generate(nEvt)};
401 } else {
402 // Make empty dataset
403 ownedGenSample = std::make_unique<RooDataSet>("emptySample","emptySample",_dependents);
404 }
405 }
406
407 _genSample = ownedGenSample.get();
408
409 //} else if (asciiFilePat && &asciiFilePat) { //warning: the address of 'asciiFilePat' will always evaluate as 'true'
410 } else if (asciiFilePat) {
411
412 // Load sample from ASCII file
413 char asciiFile[1024] ;
414 snprintf(asciiFile,1024,asciiFilePat,nSamples) ;
415 RooArgList depList(_allDependents) ;
416 ownedGenSample = std::unique_ptr<RooDataSet>{RooDataSet::read(asciiFile,depList,"q")};
417 _genSample = ownedGenSample.get();
418
419 } else {
420
421 // Load sample from internal list
422 _genSample = static_cast<RooDataSet*>(_genDataList.At(nSamples)) ;
423 existingData = true ;
424 if (!_genSample) {
425 oocoutW(_fitModel,Generation) << "RooMCStudy::run: WARNING: Sample #" << nSamples << " not loaded, skipping" << std::endl ;
426 continue ;
427 }
428 }
429
430 // Save number of generated events
431 _ngenVar->setVal(_genSample->sumEntries()) ;
432
433 // Call module between generation and fitting hook
434 for (RooAbsMCStudyModule *mod : _modList) {
435 mod->processBetweenGenAndFit(nSamples) ;
436 }
437
438 if (DoFit) fitSample(_genSample) ;
439
440 // Call module between generation and fitting hook
441 for (RooAbsMCStudyModule *mod : _modList) {
442 mod->processAfterFit(nSamples) ;
443 }
444
445 // Optionally write to ascii file
446 if (doGenerate && asciiFilePat && *asciiFilePat) {
447 char asciiFile[1024] ;
448 snprintf(asciiFile,1024,asciiFilePat,nSamples) ;
449 if (RooDataSet* unbinnedData = dynamic_cast<RooDataSet*>(_genSample)) {
450 unbinnedData->write(asciiFile) ;
451 } else {
452 coutE(InputArguments) << "RooMCStudy::run(" << GetName() << ") ERROR: ASCII writing of binned datasets is not supported" << std::endl ;
453 }
454 }
455
456 // Add to list or delete
457 if (!existingData) {
458 if (keepGenData) {
459 _genDataList.Add(ownedGenSample.release()) ;
460 }
461 }
462 }
463
464 for (RooAbsMCStudyModule *mod : _modList) {
465 if (RooDataSet* auxData = mod->finalizeRun()) {
466 _fitParData->merge(auxData) ;
467 }
468 }
469
470 _canAddFitResults = false ;
471
472 if (_genParData) {
473 for(RooAbsArg * arg : *_genParData->get()) {
474 _genParData->changeObservableName(arg->GetName(),(std::string(arg->GetName()) + "_gen").c_str());
475 }
476
477 _fitParData->merge(_genParData.get());
478 }
479
480 if (DoFit) calcPulls() ;
481
482 if (_silence) {
484 }
485
486 return false ;
487}
488
489
490
491
492
493
494////////////////////////////////////////////////////////////////////////////////
495/// Generate and fit 'nSamples' samples of 'nEvtPerSample' events.
496/// If keepGenData is set, all generated data sets will be kept in memory and can be accessed
497/// later via genData().
498///
499/// Data sets will be written out in ascii form if the pattern string is supplied.
500/// The pattern, which is a template for snprintf, should look something like "data/toymc_%04d.dat"
501/// and should contain one integer field that encodes the sample serial number.
502///
503
504bool RooMCStudy::generateAndFit(Int_t nSamples, Int_t nEvtPerSample, bool keepGenData, const char* asciiFilePat)
505{
506 // Clear any previous data in memory
507 _fitResList.Delete() ; // even though the fit results are owned by gROOT, we still want to scratch them here.
509 _fitParData->reset() ;
510
511 return run(true,true,nSamples,nEvtPerSample,keepGenData,asciiFilePat) ;
512}
513
514
515
516////////////////////////////////////////////////////////////////////////////////
517/// Generate 'nSamples' samples of 'nEvtPerSample' events.
518/// If keepGenData is set, all generated data sets will be kept in memory
519/// and can be accessed later via genData().
520///
521/// Data sets will be written out in ascii form if the pattern string is supplied.
522/// The pattern, which is a template for snprintf, should look something like "data/toymc_%04d.dat"
523/// and should contain one integer field that encodes the sample serial number.
524///
525
526bool RooMCStudy::generate(Int_t nSamples, Int_t nEvtPerSample, bool keepGenData, const char* asciiFilePat)
527{
528 // Clear any previous data in memory
530
531 return run(true,false,nSamples,nEvtPerSample,keepGenData,asciiFilePat) ;
532}
533
534
535
536////////////////////////////////////////////////////////////////////////////////
537/// Fit 'nSamples' datasets, which are read from ASCII files.
538///
539/// The ascii file pattern, which is a template for snprintf, should look something like "data/toymc_%04d.dat"
540/// and should contain one integer field that encodes the sample serial number.
541///
542
543bool RooMCStudy::fit(Int_t nSamples, const char* asciiFilePat)
544{
545 // Clear any previous data in memory
546 _fitResList.Delete() ; // even though the fit results are owned by gROOT, we still want to scratch them here.
547 _fitParData->reset() ;
548
549 return run(false,true,nSamples,0,false,asciiFilePat) ;
550}
551
552
553
554////////////////////////////////////////////////////////////////////////////////
555/// Fit 'nSamples' datasets, as supplied in 'dataSetList'
556///
557
558bool RooMCStudy::fit(Int_t nSamples, TList& dataSetList)
559{
560 // Clear any previous data in memory
561 _fitResList.Delete() ; // even though the fit results are owned by gROOT, we still want to scratch them here.
563 _fitParData->reset() ;
564
565 // Load list of data sets
566 for(auto * gset : static_range_cast<RooAbsData*>(dataSetList)) {
567 _genDataList.Add(gset) ;
568 }
569
570 return run(false,true,nSamples,0,true,nullptr) ;
571}
572
573
574
575////////////////////////////////////////////////////////////////////////////////
576/// Reset all fit parameters to the initial model
577/// parameters at the time of the RooMCStudy constructor
578
580{
582}
583
584
585
586////////////////////////////////////////////////////////////////////////////////
587/// Internal function. Performs actual fit according to specifications
588
590{
591 // Optionally bin dataset before fitting
592 std::unique_ptr<RooDataHist> ownedDataHist;
594 if (_binGenData) {
595 RooArgSet depList;
596 _fitModel->getObservables(genSample->get(), depList);
597 ownedDataHist = std::make_unique<RooDataHist>(genSample->GetName(),genSample->GetTitle(),depList,*genSample) ;
598 data = ownedDataHist.get();
599 } else {
600 data = genSample ;
601 }
602
603 RooCmdArg save = RooFit::Save() ;
605 RooCmdArg plevel = RooFit::PrintLevel(_silence ? -1 : 1) ;
606
607 RooLinkedList fitOptList(_fitOptList) ;
608 fitOptList.Add(&save) ;
609 if (!_projDeps.empty()) {
610 fitOptList.Add(&condo) ;
611 }
612 fitOptList.Add(&plevel) ;
613 return _fitModel->fitTo(*data,fitOptList);
614}
615
616
617
618////////////////////////////////////////////////////////////////////////////////
619/// Redo fit on 'current' toy sample, or if genSample is not nullptr
620/// do fit on given sample instead
621
623{
624 if (!genSample) {
625 genSample = _genSample ;
626 }
627
628 std::unique_ptr<RooFitResult> fr;
629 if (genSample->sumEntries()>0) {
630 fr = std::unique_ptr<RooFitResult>{doFit(genSample)};
631 }
632
633 return RooFit::makeOwningPtr(std::move(fr));
634}
635
636
637
638////////////////////////////////////////////////////////////////////////////////
639/// Internal method. Fit given dataset with fit model. If fit
640/// converges (TMinuit status code zero) The fit results are appended
641/// to the fit results dataset
642///
643/// If the fit option "r" is supplied, the RooFitResult
644/// objects will always be saved, regardless of the
645/// fit status. RooFitResults objects can be retrieved
646/// later via fitResult().
647///
648
650{
651 // Reset all fit parameters to their initial values
653
654 // Perform actual fit
655 bool ok ;
656 std::unique_ptr<RooFitResult> fr;
657 if (genSample->sumEntries()>0) {
658 fr = std::unique_ptr<RooFitResult>{doFit(genSample)};
659 ok = (fr->status()==0) ;
660 } else {
661 ok = false ;
662 }
663
664 // If fit converged, store parameters and NLL
665 if (ok) {
666 _nllVar->setVal(fr->minNll()) ;
667 RooArgSet tmp(_fitParams) ;
668 tmp.add(*_nllVar) ;
669 tmp.add(*_ngenVar) ;
670
671 _fitParData->add(tmp) ;
672 }
673
674 // Store fit result if requested by user
675 if (_fitOptList.FindObject("Save")) {
676 _fitResList.Add(fr.release()) ;
677 }
678
679 return !ok ;
680}
681
682
683
684////////////////////////////////////////////////////////////////////////////////
685/// Utility function to add fit result from external fit to this RooMCStudy
686/// and process its results through the standard RooMCStudy statistics gathering tools.
687/// This function allows users to run the toy MC generation and/or fitting
688/// in a distributed way and to collect and analyze the results in a RooMCStudy
689/// as if they were run locally.
690///
691/// This method is only functional if this RooMCStudy object is cleanm, i.e. it was not used
692/// to generate and/or fit any samples.
693
695{
696 if (!_canAddFitResults) {
697 oocoutE(_fitModel,InputArguments) << "RooMCStudy::addFitResult: ERROR cannot add fit results in current state" << std::endl ;
698 return true ;
699 }
700
701 // Transfer contents of fit result to fitParams ;
703
704 // If fit converged, store parameters and NLL
705 bool ok = (fr.status()==0) ;
706 if (ok) {
707 _nllVar->setVal(fr.minNll()) ;
708 RooArgSet tmp(_fitParams) ;
709 tmp.add(*_nllVar) ;
710 tmp.add(*_ngenVar) ;
711 _fitParData->add(tmp) ;
712 }
713
714 // Store fit result if requested by user
715 if (_fitOptList.FindObject("Save")) {
716 _fitResList.Add((TObject*)&fr) ;
717 }
718
719 return false ;
720}
721
722
723
724////////////////////////////////////////////////////////////////////////////////
725/// Calculate the pulls for all fit parameters in
726/// the fit results data set, and add them to that dataset.
727
729{
730 for (auto it = _fitParams.begin(); it != _fitParams.end(); ++it) {
731 const auto par = static_cast<RooRealVar*>(*it);
732 _fitParData->addColumn(*std::unique_ptr<RooErrorVar>{par->errorVar()});
733
734 TString name(par->GetName());
735 TString title(par->GetTitle());
736 name.Append("pull") ;
737 title.Append(" Pull") ;
738
739 if (!par->hasError(false)) {
740 coutW(Generation) << "Fit parameter '" << par->GetName() << "' does not have an error."
741 " A pull distribution cannot be generated. This might be caused by the parameter being constant or"
742 " because the fits were not run." << std::endl;
743 continue;
744 }
745
746 // First look in fitParDataset to see if per-experiment generated value has been stored
747 auto genParOrig = static_cast<RooAbsReal*>(_fitParData->get()->find(Form("%s_gen",par->GetName())));
748 if (genParOrig && _perExptGenParams) {
749
750 RooPullVar pull(name,title,*par,*genParOrig) ;
751 _fitParData->addColumn(pull,false) ;
752
753 } else {
754 // If not use fixed generator value
755 genParOrig = static_cast<RooAbsReal*>(_genInitParams.find(par->GetName()));
756
757 if (!genParOrig) {
758 std::size_t index = it - _fitParams.begin();
759 genParOrig = index < _genInitParams.size() ?
760 static_cast<RooAbsReal*>(_genInitParams[index]) :
761 nullptr;
762
763 if (genParOrig) {
764 coutW(Generation) << "The fit parameter '" << par->GetName() << "' is not in the model that was used to generate toy data. "
765 "The parameter '" << genParOrig->GetName() << "'=" << genParOrig->getVal() << " was found at the same position in the generator model."
766 " It will be used to compute pulls."
767 "\nIf this is not desired, the parameters of the generator model need to be renamed or reordered." << std::endl;
768 }
769 }
770
771 if (genParOrig) {
772 std::unique_ptr<RooAbsReal> genPar(static_cast<RooAbsReal*>(genParOrig->Clone("truth")));
773 RooPullVar pull(name,title,*par,*genPar);
774
775 _fitParData->addColumn(pull,false) ;
776 } else {
777 coutE(Generation) << "Cannot generate pull distribution for the fit parameter '" << par->GetName() << "'."
778 "\nNo similar parameter was found in the set of parameters that were used to generate toy data." << std::endl;
779 }
780 }
781 }
782}
783
784
785
786
787////////////////////////////////////////////////////////////////////////////////
788/// Return a RooDataSet containing the post-fit parameters of each toy cycle.
789/// This dataset also contains any additional output that was generated
790/// by study modules that were added to this RooMCStudy.
791/// By default, the two following variables are added (apart from fit parameters):
792/// - NLL: The value of the negative log-likelihood for each run.
793/// - ngen: Number of events generated for each run.
795{
796 if (_canAddFitResults) {
797 calcPulls() ;
798 _canAddFitResults = false ;
799 }
800
801 return *_fitParData ;
802}
803
804
805
806////////////////////////////////////////////////////////////////////////////////
807/// Return an argset with the fit parameters for the given sample number
808///
809/// NB: The fit parameters are only stored for successful fits,
810/// thus the maximum sampleNum can be less that the number
811/// of generated samples and if so, the indices will
812/// be out of synch with genData() and fitResult()
813
814const RooArgSet* RooMCStudy::fitParams(Int_t sampleNum) const
815{
816 // Check if sampleNum is in range
817 if (sampleNum<0 || sampleNum>=_fitParData->numEntries()) {
818 oocoutE(_fitModel,InputArguments) << "RooMCStudy::fitParams: ERROR, invalid sample number: " << sampleNum << std::endl ;
819 return nullptr ;
820 }
821
822 return _fitParData->get(sampleNum) ;
823}
824
825
826
827////////////////////////////////////////////////////////////////////////////////
828/// Return the RooFitResult of the fit with the given run number.
829///
830/// \note Fit results are not saved by default. This requires passing `FitOptions(Save(), ...)`
831/// to the constructor.
833{
834 // Check if sampleNum is in range
835 if (sampleNum<0 || sampleNum>=_fitResList.GetSize()) {
836 oocoutE(_fitModel,InputArguments) << "RooMCStudy::fitResult: ERROR, invalid sample number: " << sampleNum << std::endl ;
837 return nullptr ;
838 }
839
840 // Retrieve fit result object
841 const RooFitResult* fr = static_cast<RooFitResult*>(_fitResList.At(sampleNum)) ;
842 if (fr) {
843 return fr ;
844 } else {
845 oocoutE(_fitModel,InputArguments) << "RooMCStudy::fitResult: ERROR, no fit result saved for sample "
846 << sampleNum << ", did you use the 'r; fit option?" << std::endl ;
847 }
848 return nullptr ;
849}
850
851
852
853////////////////////////////////////////////////////////////////////////////////
854/// Return the given generated dataset. This method will only return datasets
855/// if during the run cycle it was indicated that generator data should be saved.
856
858{
859 // Check that generated data was saved
860 if (_genDataList.GetSize()==0) {
861 oocoutE(_fitModel,InputArguments) << "RooMCStudy::genData() ERROR, generated data was not saved" << std::endl ;
862 return nullptr ;
863 }
864
865 // Check if sampleNum is in range
866 if (sampleNum<0 || sampleNum>=_genDataList.GetSize()) {
867 oocoutE(_fitModel,InputArguments) << "RooMCStudy::genData() ERROR, invalid sample number: " << sampleNum << std::endl ;
868 return nullptr ;
869 }
870
871 return static_cast<RooAbsData*>(_genDataList.At(sampleNum)) ;
872}
873
874
875
876////////////////////////////////////////////////////////////////////////////////
877/// Plot the distribution of fitted values of a parameter. The parameter shown is the one from which the RooPlot
878/// was created, e.g.
879///
880/// RooPlot* frame = param.frame(100,-10,10) ;
881/// mcstudy.paramOn(frame,LineStyle(kDashed)) ;
882///
883/// Any named arguments passed to plotParamOn() are forwarded to the underlying plotOn() call
884
885RooPlot* RooMCStudy::plotParamOn(RooPlot* frame, const RooCmdArg& arg1, const RooCmdArg& arg2, const RooCmdArg& arg3, const RooCmdArg& arg4,
886 const RooCmdArg& arg5, const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8)
887{
888 _fitParData->plotOn(frame,arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8) ;
889 return frame ;
890}
891
892
893
894////////////////////////////////////////////////////////////////////////////////
895/// Plot the distribution of the fitted value of the given parameter on a newly created frame.
896///
897/// <table>
898/// <tr><th> Optional arguments <th>
899/// <tr><td> FrameRange(double lo, double hi) <td> Set range of frame to given specification
900/// <tr><td> FrameBins(int bins) <td> Set default number of bins of frame to given number
901/// <tr><td> Frame() <td> Pass supplied named arguments to RooAbsRealLValue::frame() function. See there
902/// for list of allowed arguments
903/// </table>
904/// If no frame specifications are given, the AutoRange() feature will be used to set the range
905/// Any other named argument is passed to the RooAbsData::plotOn() call. See that function for allowed options
906
907RooPlot* RooMCStudy::plotParam(const char* paramName, const RooCmdArg& arg1, const RooCmdArg& arg2, const RooCmdArg& arg3, const RooCmdArg& arg4,
908 const RooCmdArg& arg5, const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8)
909{
910
911 // Find parameter in fitParDataSet
912 RooRealVar* param = static_cast<RooRealVar*>(_fitParData->get()->find(paramName)) ;
913 if (!param) {
914 oocoutE(_fitModel,InputArguments) << "RooMCStudy::plotParam: ERROR: no parameter defined with name " << paramName << std::endl ;
915 return nullptr ;
916 }
917
918 // Forward to implementation below
919 return plotParam(*param,arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8) ;
920}
921
922
923
924////////////////////////////////////////////////////////////////////////////////
925/// Plot the distribution of the fitted value of the given parameter on a newly created frame.
926/// \copydetails RooMCStudy::plotParam(const char* paramName, const RooCmdArg& arg1, const RooCmdArg& arg2, const RooCmdArg& arg3, const RooCmdArg& arg4,
927/// const RooCmdArg& arg5, const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8)
928
929RooPlot* RooMCStudy::plotParam(const RooRealVar& param, const RooCmdArg& arg1, const RooCmdArg& arg2, const RooCmdArg& arg3, const RooCmdArg& arg4,
930 const RooCmdArg& arg5, const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8)
931{
932 // Stuff all arguments in a list
933 RooLinkedList cmdList;
934 cmdList.Add(const_cast<RooCmdArg*>(&arg1)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg2)) ;
935 cmdList.Add(const_cast<RooCmdArg*>(&arg3)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg4)) ;
936 cmdList.Add(const_cast<RooCmdArg*>(&arg5)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg6)) ;
937 cmdList.Add(const_cast<RooCmdArg*>(&arg7)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg8)) ;
938
939 RooPlot* frame = makeFrameAndPlotCmd(param, cmdList) ;
940 if (frame) {
941 _fitParData->plotOn(frame, cmdList) ;
942 }
943
944 return frame ;
945}
946
947
948
949////////////////////////////////////////////////////////////////////////////////
950/// Plot the distribution of the -log(L) values on a newly created frame.
951///
952/// <table>
953/// <tr><th> Optional arguments <th>
954/// <tr><td> FrameRange(double lo, double hi) <td> Set range of frame to given specification
955/// <tr><td> FrameBins(int bins) <td> Set default number of bins of frame to given number
956/// <tr><td> Frame() <td> Pass supplied named arguments to RooAbsRealLValue::frame() function. See there
957/// for list of allowed arguments
958/// </table>
959///
960/// If no frame specifications are given, the AutoRange() feature will be used to set the range.
961/// Any other named argument is passed to the RooAbsData::plotOn() call. See that function for allowed options
962
964 const RooCmdArg& arg3, const RooCmdArg& arg4,
965 const RooCmdArg& arg5, const RooCmdArg& arg6,
966 const RooCmdArg& arg7, const RooCmdArg& arg8)
967{
968 return plotParam(*_nllVar,arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8) ;
969}
970
971
972
973////////////////////////////////////////////////////////////////////////////////
974/// Plot the distribution of the fit errors for the specified parameter on a newly created frame.
975///
976/// <table>
977/// <tr><th> Optional arguments <th>
978/// <tr><td> FrameRange(double lo, double hi) <td> Set range of frame to given specification
979/// <tr><td> FrameBins(int bins) <td> Set default number of bins of frame to given number
980/// <tr><td> Frame() <td> Pass supplied named arguments to RooAbsRealLValue::frame() function. See there
981/// for list of allowed arguments
982/// </table>
983///
984/// If no frame specifications are given, the AutoRange() feature will be used to set a default range.
985/// Any other named argument is passed to the RooAbsData::plotOn() call. See that function for allowed options.
986
987RooPlot* RooMCStudy::plotError(const RooRealVar& param, const RooCmdArg& arg1, const RooCmdArg& arg2,
988 const RooCmdArg& arg3, const RooCmdArg& arg4,
989 const RooCmdArg& arg5, const RooCmdArg& arg6,
990 const RooCmdArg& arg7, const RooCmdArg& arg8)
991{
992 if (_canAddFitResults) {
993 calcPulls() ;
994 _canAddFitResults=false ;
995 }
996
997 std::unique_ptr<RooErrorVar> evar{param.errorVar()};
998 std::unique_ptr<RooAbsArg> evar_rrv{evar->createFundamental()};
999 RooPlot* frame = plotParam(static_cast<RooRealVar&>(*evar_rrv),arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8) ;
1000
1001 // To make sure the frame has no dangling pointer to evar_rrv.
1003
1004 return frame ;
1005}
1006
1007namespace {
1008
1009// Fits a Gaussian to the pull distribution, plots the fit and prints the fit
1010// parameters on the canvas. Implementation detail of RooMCStudy::plotPull().
1011void fitGaussToPulls(RooPlot& frame, RooDataSet& fitParData)
1012{
1013 // Build the Gaussian fit mode for the pulls, then fit it and plot it. We
1014 // have to use the RooWorkspace factory here, because different from the
1015 // RooMCStudy class, the RooGaussian is not in RooFitCore.
1016 RooWorkspace ws;
1017 auto plotVar = frame.getPlotVar();
1018 const std::string plotVarName = plotVar->GetName();
1019 ws.import(*plotVar);
1020 ws.factory("Gaussian::pullGauss(" + plotVarName + ", pullMean[0.0, -10.0, 10.0], pullSigma[1.0, 0.1, 5.0])");
1021
1022 RooRealVar& pullMean = *ws.var("pullMean");
1023 RooRealVar& pullSigma = *ws.var("pullSigma");
1024 RooAbsPdf& pullGauss = *ws.pdf("pullGauss");
1025
1026 pullGauss.fitTo(fitParData, RooFit::Minos(false), RooFit::PrintLevel(-1)) ;
1027 pullGauss.plotOn(&frame) ;
1028
1029 // Instead of using paramOn() without command arguments to plot the fit
1030 // parameters, we are building the parameter label ourselves for more
1031 // flexibility and pass this together with an appropriate layout
1032 // parametrization to paramOn().
1033 const int sigDigits = 2;
1034 const char * options = "ELU";
1035 std::stringstream ss;
1036 ss << "Fit parameters:\n"
1037 << "#mu: " << *std::unique_ptr<TString>{pullMean.format(sigDigits, options)}
1038 << "\n#sigma: " << *std::unique_ptr<TString>{pullSigma.format(sigDigits, options)};
1039 // We set the parameters constant to disable the default label. Still, we
1040 // use param() on as a wrapper for the text box generation.
1041 pullMean.setConstant(true);
1042 pullSigma.setConstant(true);
1043 pullGauss.paramOn(&frame, RooFit::Label(ss.str().c_str()), RooFit::Layout(0.60, 0.9, 0.9));
1044}
1045
1046} // namespace
1047
1048
1049////////////////////////////////////////////////////////////////////////////////
1050/// Plot the distribution of pull values for the specified parameter on a newly created frame. If asymmetric
1051/// errors are calculated in the fit (by MINOS) those will be used in the pull calculation.
1052///
1053/// If the parameters of the models for generation and fit differ, simple heuristics are used to find the
1054/// corresponding parameters:
1055/// - Parameters have the same name: They will be used to compute pulls.
1056/// - Parameters have different names: The position of the fit parameter in the set of fit parameters will be
1057/// computed. The parameter at the same position in the set of generator parameters will be used.
1058///
1059/// Further options:
1060/// <table>
1061/// <tr><th> Arguments <th> Effect
1062/// <tr><td> FrameRange(double lo, double hi) <td> Set range of frame to given specification
1063/// <tr><td> FrameBins(int bins) <td> Set default number of bins of frame to given number
1064/// <tr><td> Frame() <td> Pass supplied named arguments to RooAbsRealLValue::frame() function. See there
1065/// for list of allowed arguments
1066/// <tr><td> FitGauss(bool flag) <td> Add a gaussian fit to the frame
1067/// </table>
1068///
1069/// If no frame specifications are given, the AutoSymRange() feature will be used to set a default range.
1070/// Any other named argument is passed to the RooAbsData::plotOn(). See that function for allowed options.
1071///
1072/// If you want to have more control over the Gaussian fit to the pull
1073/// distribution, you can also do it after the call to plotPull():
1074///
1075/// ~~~ {.cpp}
1076/// RooPlot *frame = mcstudy->plotPull(myVariable, RooFit::Bins(40), RooFit::FitGauss(false));
1077/// RooRealVar pullMean("pullMean","Mean of pull",0,-10,10) ;
1078/// RooRealVar pullSigma("pullSigma","Width of pull",1,0.1,5) ;
1079/// pullMean.setPlotLabel("pull #mu"); // optional (to get nicer plot labels if you want)
1080/// pullSigma.setPlotLabel("pull #sigma"); // optional
1081/// RooGaussian pullGauss("pullGauss","Gaussian of pull", *frame->getPlotVar(), pullMean, pullSigma);
1082/// pullGauss.fitTo(const_cast<RooDataSet&>(mcstudy->fitParDataSet()),
1083/// RooFit::Minos(0), RooFit::PrintLevel(-1)) ;
1084/// pullGauss.plotOn(frame) ;
1085/// pullGauss.paramOn(frame, RooFit::Layout(0.65, 0.9, 0.9)); // optionally specify label position (xmin, xmax, ymax)
1086/// ~~~
1087
1088RooPlot* RooMCStudy::plotPull(const RooRealVar& param, const RooCmdArg& arg1, const RooCmdArg& arg2,
1089 const RooCmdArg& arg3, const RooCmdArg& arg4,
1090 const RooCmdArg& arg5, const RooCmdArg& arg6,
1091 const RooCmdArg& arg7, const RooCmdArg& arg8)
1092{
1093 // Stuff all arguments in a list
1094 RooLinkedList cmdList;
1095 cmdList.Add(const_cast<RooCmdArg*>(&arg1)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg2)) ;
1096 cmdList.Add(const_cast<RooCmdArg*>(&arg3)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg4)) ;
1097 cmdList.Add(const_cast<RooCmdArg*>(&arg5)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg6)) ;
1098 cmdList.Add(const_cast<RooCmdArg*>(&arg7)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg8)) ;
1099
1100 TString name(param.GetName());
1101 TString title(param.GetTitle());
1102 name.Append("pull") ; title.Append(" Pull") ;
1103 RooRealVar pvar(name,title,-100,100) ;
1104 pvar.setBins(100) ;
1105
1106
1107 RooPlot* frame = makeFrameAndPlotCmd(pvar, cmdList, true) ;
1108 if (frame) {
1109
1110 // Pick up optional FitGauss command from list
1111 RooCmdConfig pc("RooMCStudy::plotPull(" + std::string(_genModel->GetName()) + ")");
1112 pc.defineInt("fitGauss","FitGauss",0,0) ;
1113 pc.allowUndefined() ;
1114 pc.process(cmdList) ;
1115 bool fitGauss=pc.getInt("fitGauss") ;
1116
1117 // Pass stripped command list to plotOn()
1118 RooCmdConfig::stripCmdList(cmdList,"FitGauss") ;
1119 const bool success = _fitParData->plotOn(frame,cmdList) ;
1120
1121 if (!success) {
1122 coutF(Plotting) << "No pull distribution for the parameter '" << param.GetName() << "'. Check logs for errors." << std::endl;
1123 return frame;
1124 }
1125
1126 // Add Gaussian fit if requested
1127 if (fitGauss) {
1128 fitGaussToPulls(*frame, *_fitParData);
1129 }
1130
1131 // To make sure the frame has no dangling pointer to pvar.
1133 }
1134 return frame;
1135}
1136
1137
1138
1139////////////////////////////////////////////////////////////////////////////////
1140/// Internal function. Construct RooPlot from given parameter and modify the list of named
1141/// arguments 'cmdList' to only contain the plot arguments that should be forwarded to
1142/// RooAbsData::plotOn()
1143
1144RooPlot* RooMCStudy::makeFrameAndPlotCmd(const RooRealVar& param, RooLinkedList& cmdList, bool symRange) const
1145{
1146 // Select the frame-specific commands
1147 RooCmdConfig pc("RooMCStudy::plotParam(" + std::string(_genModel->GetName()) + ")");
1148 pc.defineInt("nbins","Bins",0,0) ;
1149 pc.defineDouble("xlo","Range",0,0) ;
1150 pc.defineDouble("xhi","Range",1,0) ;
1151 pc.defineInt("dummy","FrameArgs",0,0) ;
1152 pc.defineMutex("Bins","FrameArgs") ;
1153 pc.defineMutex("Range","FrameArgs") ;
1154
1155 // Process and check varargs
1156 pc.allowUndefined() ;
1157 pc.process(cmdList) ;
1158 if (!pc.ok(true)) {
1159 return nullptr ;
1160 }
1161
1162 // Make frame according to specs
1163 Int_t nbins = pc.getInt("nbins") ;
1164 double xlo = pc.getDouble("xlo") ;
1165 double xhi = pc.getDouble("xhi") ;
1166 RooPlot* frame ;
1167
1168 if (pc.hasProcessed("FrameArgs")) {
1169 // Explicit frame arguments are given, pass them on
1170 RooCmdArg* frameArg = static_cast<RooCmdArg*>(cmdList.FindObject("FrameArgs")) ;
1171 frame = param.frame(frameArg->subArgs()) ;
1172 } else {
1173 // FrameBins, FrameRange or none are given, build custom frame command list
1174 RooCmdArg bins = RooFit::Bins(nbins) ;
1175 RooCmdArg range = RooFit::Range(xlo,xhi) ;
1176 RooCmdArg autoRange = symRange ? RooFit::AutoSymRange(*_fitParData,0.2) : RooFit::AutoRange(*_fitParData,0.2) ;
1177 RooLinkedList frameCmdList ;
1178
1179 if (pc.hasProcessed("Bins")) frameCmdList.Add(&bins) ;
1180 if (pc.hasProcessed("Range")) {
1181 frameCmdList.Add(&range) ;
1182 } else {
1183 frameCmdList.Add(&autoRange) ;
1184 }
1185 frame = param.frame(frameCmdList) ;
1186 }
1187
1188 // Filter frame command from list and pass on to plotOn()
1189 RooCmdConfig::stripCmdList(cmdList,"FrameArgs,Bins,Range") ;
1190
1191 return frame ;
1192}
1193
1194
1195
1196////////////////////////////////////////////////////////////////////////////////
1197/// Create a RooPlot of the -log(L) distribution in the range lo-hi
1198/// with 'nBins' bins
1199
1200RooPlot* RooMCStudy::plotNLL(double lo, double hi, Int_t nBins)
1201{
1202 RooPlot* frame = _nllVar->frame(lo,hi,nBins) ;
1203
1204 _fitParData->plotOn(frame) ;
1205 return frame ;
1206}
1207
1208
1209
1210////////////////////////////////////////////////////////////////////////////////
1211/// Create a RooPlot of the distribution of the fitted errors of the given parameter.
1212/// The frame is created with a range [lo,hi] and plotted data will be binned in 'nbins' bins
1213
1214RooPlot* RooMCStudy::plotError(const RooRealVar& param, double lo, double hi, Int_t nbins)
1215{
1216 if (_canAddFitResults) {
1217 calcPulls() ;
1218 _canAddFitResults=false ;
1219 }
1220
1221 std::unique_ptr<RooErrorVar> evar{param.errorVar()};
1222 RooPlot* frame = evar->frame(lo,hi,nbins) ;
1223 _fitParData->plotOn(frame) ;
1224
1225 return frame ;
1226}
1227
1228
1229
1230////////////////////////////////////////////////////////////////////////////////
1231/// Create a RooPlot of the pull distribution for the given
1232/// parameter. The range lo-hi is plotted in nbins. If fitGauss is
1233/// set, an unbinned ML fit of the distribution to a Gaussian p.d.f
1234/// is performed. The fit result is overlaid on the returned RooPlot
1235/// and a box with the fitted mean and sigma is added.
1236///
1237/// If the parameters of the models for generation and fit differ, simple heuristics are used to find the
1238/// corresponding parameters:
1239/// - Parameters have the same name: They will be used to compute pulls.
1240/// - Parameters have different names: The position of the fit parameter in the set of fit parameters will be
1241/// computed. The parameter at the same position in the set of generator parameters will be used.
1242
1243RooPlot* RooMCStudy::plotPull(const RooRealVar& param, double lo, double hi, Int_t nbins, bool fitGauss)
1244{
1245 if (_canAddFitResults) {
1246 calcPulls() ;
1247 _canAddFitResults=false ;
1248 }
1249
1250 TString name(param.GetName());
1251 TString title(param.GetTitle());
1252 name.Append("pull") ; title.Append(" Pull") ;
1253 RooRealVar pvar(name,title,lo,hi) ;
1254 pvar.setBins(nbins) ;
1255
1256 RooPlot* frame = pvar.frame() ;
1257 const bool success = _fitParData->plotOn(frame);
1258
1259 if (!success) {
1260 coutF(Plotting) << "No pull distribution for the parameter '" << param.GetName() << "'. Check logs for errors." << std::endl;
1261 return frame;
1262 }
1263
1264 if (fitGauss) {
1265 fitGaussToPulls(*frame, *_fitParData);
1266 }
1267
1268 return frame ;
1269}
1270
1271
1272////////////////////////////////////////////////////////////////////////////////
1273/// If one of the TObject we have a referenced to is deleted, remove the
1274/// reference.
1275
1277{
1281 if (_ngenVar.get() == obj) _ngenVar.reset();
1282
1283 if (_fitParData) _fitParData->RecursiveRemove(obj);
1284 if (_fitParData.get() == obj) _fitParData.reset();
1285
1286 if (_genParData) _genParData->RecursiveRemove(obj);
1287 if (_genParData.get() == obj) _genParData.reset();
1288}
1289
#define coutI(a)
#define oocoutW(o, a)
#define coutW(a)
#define coutF(a)
#define oocoutE(o, a)
#define oocoutI(o, a)
#define coutE(a)
#define ooccoutI(o, a)
#define ooccoutP(o, a)
#define oocoutP(o, a)
#define ClassImp(name)
Definition Rtypes.h:377
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t index
char name[80]
Definition TGX11.cxx:110
#define hi
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition TString.cxx:2489
#define snprintf
Definition civetweb.c:1540
Common abstract base class for objects that represent a value and a "shape" in RooFit.
Definition RooAbsArg.h:77
RooFit::OwningPtr< RooArgSet > getParameters(const RooAbsData *data, bool stripDisconnected=true) const
Create a list of leaf nodes in the arg tree starting with ourself as top node that don't match any of...
RooFit::OwningPtr< RooArgSet > getObservables(const RooArgSet &set, bool valueOnly=true) const
Given a set of possible observables, return the observables that this PDF depends on.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
void setAttribAll(const Text_t *name, bool value=true)
Set given attribute in each element of the collection by calling each elements setAttribute() functio...
void assign(const RooAbsCollection &other) const
Sets the value, cache and constant attribute of any argument in our set that also appears in the othe...
Storage_t::size_type size() const
TIterator begin()
TIterator end() and range-based for loops.")
RooAbsArg * find(const char *name) const
Find object with given name in list.
Abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:57
virtual double sumEntries() const =0
Return effective number of entries in dataset, i.e., sum all weights.
virtual const RooArgSet * get() const
Definition RooAbsData.h:101
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Base class for add-on modules to RooMCStudy that can perform additional calculations on each generate...
bool doInitializeInstance(RooMCStudy &)
Initializer method called upon attachment to given RooMCStudy object.
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
virtual double expectedEvents(const RooArgSet *nset) const
Return expected number of events to be used in calculation of extended likelihood.
Int_t * randomizeProtoOrder(Int_t nProto, Int_t nGen, bool resample=false) const
Return lookup table with randomized order for nProto prototype events.
RooPlot * plotOn(RooPlot *frame, const RooCmdArg &arg1={}, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}, const RooCmdArg &arg9={}, const RooCmdArg &arg10={}) const override
Helper calling plotOn(RooPlot*, RooLinkedList&) const.
Definition RooAbsPdf.h:124
virtual RooPlot * paramOn(RooPlot *frame, const RooCmdArg &arg1={}, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={})
Add a box with parameter values (and errors) to the specified frame.
RooFit::OwningPtr< RooFitResult > fitTo(RooAbsData &data, CmdArgs_t const &... cmdArgs)
Fit PDF to given dataset.
Definition RooAbsPdf.h:157
RooArgSet * getAllConstraints(const RooArgSet &observables, RooArgSet &constrainedParams, bool stripDisconnected=true, bool removeConstraintsFromPdf=false) const
This helper function finds and collects all constraints terms of all component p.d....
virtual RooFit::OwningPtr< RooDataHist > generateBinned(const RooArgSet &whatVars, double nEvents, const RooCmdArg &arg1, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}) const
As RooAbsPdf::generateBinned(const RooArgSet&, const RooCmdArg&,const RooCmdArg&, const RooCmdArg&,...
Definition RooAbsPdf.h:110
virtual RooAbsGenContext * genContext(const RooArgSet &vars, const RooDataSet *prototype=nullptr, const RooArgSet *auxProto=nullptr, bool verbose=false) const
Interface function to create a generator context from a p.d.f.
void setConstant(bool value=true)
RooPlot * frame(const RooCmdArg &arg1, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}) const
Create a new RooPlot on the heap with a drawing frame initialized for this object,...
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
RooFit::OwningPtr< RooAbsArg > createFundamental(const char *newname=nullptr) const override
Create a RooRealVar fundamental object with our properties.
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:22
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:55
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition RooArgSet.h:178
Named container for two doubles, two integers two object points and three string pointers that can be...
Definition RooCmdArg.h:26
RooLinkedList const & subArgs() const
Return list of sub-arguments in this RooCmdArg.
Definition RooCmdArg.h:52
TObject * Clone(const char *newName=nullptr) const override
Make a clone of an object using the Streamer facility.
Definition RooCmdArg.h:57
Configurable parser for RooCmdArg named arguments.
void defineMutex(const char *head, Args_t &&... tail)
Define arguments where any pair is mutually exclusive.
bool process(const RooCmdArg &arg)
Process given RooCmdArg.
bool hasProcessed(const char *cmdName) const
Return true if RooCmdArg with name 'cmdName' has been processed.
double getDouble(const char *name, double defaultValue=0.0) const
Return double property registered with name 'name'.
bool defineDouble(const char *name, const char *argName, int doubleNum, double defValue=0.0)
Define double property name 'name' mapped to double in slot 'doubleNum' in RooCmdArg with name argNam...
static void stripCmdList(RooLinkedList &cmdList, const char *cmdsToPurge)
Utility function that strips command names listed (comma separated) in cmdsToPurge from cmdList.
RooArgSet * getSet(const char *name, RooArgSet *set=nullptr) const
Return RooArgSet property registered with name 'name'.
bool defineSet(const char *name, const char *argName, int setNum, const RooArgSet *set=nullptr)
Define TObject property name 'name' mapped to object in slot 'setNum' in RooCmdArg with name argName ...
bool ok(bool verbose) const
Return true of parsing was successful.
bool defineObject(const char *name, const char *argName, int setNum, const TObject *obj=nullptr, bool isArray=false)
Define TObject property name 'name' mapped to object in slot 'setNum' in RooCmdArg with name argName ...
bool defineInt(const char *name, const char *argName, int intNum, int defValue=0)
Define integer property name 'name' mapped to integer in slot 'intNum' in RooCmdArg with name argName...
void allowUndefined(bool flag=true)
If flag is true the processing of unrecognized RooCmdArgs is not considered an error.
int getInt(const char *name, int defaultValue=0) const
Return integer property registered with name 'name'.
TObject * getObject(const char *name, TObject *obj=nullptr) const
Return TObject property registered with name 'name'.
Container class to hold unbinned data.
Definition RooDataSet.h:57
const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event 'index'.
static RooDataSet * read(const char *filename, const RooArgList &variables, const char *opts="", const char *commonPath="", const char *indexCatName=nullptr)
Read data from a text file and create a dataset from it.
RooFitResult is a container class to hold the input and output of a PDF fit to a dataset.
const RooArgList & floatParsFinal() const
Return list of floating parameters after fit.
Int_t status() const
Return MINUIT status code.
double minNll() const
Return minimized -log(L) value.
Collection class for internal use, storing a collection of RooAbsArg pointers in a doubly linked list...
Int_t GetSize() const
TObject * At(int index) const
Return object stored in sequential position given by index.
void RecursiveRemove(TObject *obj) override
If one of the TObject we have a referenced to is deleted, remove the reference.
void Delete(Option_t *o=nullptr) override
Remove all elements in collection and delete all elements NB: Collection does not own elements,...
virtual void Add(TObject *arg)
TObject * FindObject(const char *name) const override
Return pointer to object with given name.
Helper class to facilitate Monte Carlo studies such as 'goodness-of-fit' studies, that involve fittin...
Definition RooMCStudy.h:32
bool addFitResult(const RooFitResult &fr)
Utility function to add fit result from external fit to this RooMCStudy and process its results throu...
RooAbsData * _genSample
Currently generated sample.
Definition RooMCStudy.h:107
RooPlot * makeFrameAndPlotCmd(const RooRealVar &param, RooLinkedList &cmdList, bool symRange=false) const
Internal function.
RooArgSet _projDeps
List of projected dependents in fit.
Definition RooMCStudy.h:113
RooArgSet _genParams
List of actual generator parameters.
Definition RooMCStudy.h:111
const RooArgSet * fitParams(Int_t sampleNum) const
Return an argset with the fit parameters for the given sample number.
void calcPulls()
Calculate the pulls for all fit parameters in the fit results data set, and add them to that dataset.
~RooMCStudy() override
RooPlot * plotPull(const RooRealVar &param, const RooCmdArg &arg1, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={})
Plot the distribution of pull values for the specified parameter on a newly created frame.
RooArgSet _dependents
List of dependents.
Definition RooMCStudy.h:118
bool _verboseGen
Verbose generation?
Definition RooMCStudy.h:137
std::list< RooAbsMCStudyModule * > _modList
List of additional study modules ;.
Definition RooMCStudy.h:141
std::unique_ptr< RooDataSet > _genParData
Definition RooMCStudy.h:128
RooArgSet _genInitParams
List of original generator parameters.
Definition RooMCStudy.h:110
TList _fitResList
Definition RooMCStudy.h:127
double _nExpGen
Definition RooMCStudy.h:133
bool fitSample(RooAbsData *genSample)
Internal method.
std::unique_ptr< RooDataSet > _fitParData
Definition RooMCStudy.h:129
bool generate(Int_t nSamples, Int_t nEvtPerSample=0, bool keepGenData=false, const char *asciiFilePat=nullptr)
Generate 'nSamples' samples of 'nEvtPerSample' events.
RooPlot * plotParamOn(RooPlot *frame, const RooCmdArg &arg1={}, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={})
Plot the distribution of fitted values of a parameter.
bool _extendedGen
Definition RooMCStudy.h:131
const RooDataSet * _genProtoData
Generator prototype data set.
Definition RooMCStudy.h:112
bool _canAddFitResults
Allow adding of external fit results?
Definition RooMCStudy.h:136
const RooFitResult * fitResult(Int_t sampleNum) const
Return the RooFitResult of the fit with the given run number.
RooFit::OwningPtr< RooFitResult > doFit(RooAbsData *genSample)
Internal function. Performs actual fit according to specifications.
std::unique_ptr< RooAbsGenContext > _constrGenContext
Generator context for constraints p.d.f.
Definition RooMCStudy.h:116
bool _perExptGenParams
Do generation parameter change per event?
Definition RooMCStudy.h:138
bool _binGenData
Definition RooMCStudy.h:132
bool _silence
Silent running mode?
Definition RooMCStudy.h:139
RooArgSet _fitParams
List of actual fit parameters.
Definition RooMCStudy.h:122
std::unique_ptr< RooAbsGenContext > _genContext
Generator context.
Definition RooMCStudy.h:109
RooPlot * plotNLL(const RooCmdArg &arg1={}, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={})
Plot the distribution of the -log(L) values on a newly created frame.
RooMCStudy(const RooAbsPdf &model, const RooArgSet &observables, const RooCmdArg &arg1={}, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={})
Construct Monte Carlo Study Manager.
RooFit::OwningPtr< RooFitResult > refit(RooAbsData *genSample=nullptr)
Redo fit on 'current' toy sample, or if genSample is not nullptr do fit on given sample instead.
RooAbsData * genData(Int_t sampleNum) const
Return the given generated dataset.
void RecursiveRemove(TObject *obj) override
If one of the TObject we have a referenced to is deleted, remove the reference.
RooAbsPdf * _genModel
Generator model.
Definition RooMCStudy.h:108
const RooDataSet & fitParDataSet()
Return a RooDataSet containing the post-fit parameters of each toy cycle.
std::unique_ptr< RooRealVar > _nllVar
Definition RooMCStudy.h:123
RooLinkedList _fitOptList
Definition RooMCStudy.h:130
std::unique_ptr< RooAbsPdf > _constrPdf
Constraints p.d.f.
Definition RooMCStudy.h:115
RooArgSet _allDependents
List of generate + prototype dependents.
Definition RooMCStudy.h:119
bool run(bool generate, bool fit, Int_t nSamples, Int_t nEvtPerSample, bool keepGenData, const char *asciiFilePat)
Run engine method.
void resetFitParams()
Reset all fit parameters to the initial model parameters at the time of the RooMCStudy constructor.
RooAbsPdf * _fitModel
Fit model.
Definition RooMCStudy.h:120
bool fit(Int_t nSamples, const char *asciiFilePat)
Fit 'nSamples' datasets, which are read from ASCII files.
bool generateAndFit(Int_t nSamples, Int_t nEvtPerSample=0, bool keepGenData=false, const char *asciiFilePat=nullptr)
Generate and fit 'nSamples' samples of 'nEvtPerSample' events.
TList _genDataList
Definition RooMCStudy.h:126
bool _randProto
Definition RooMCStudy.h:134
void addModule(RooAbsMCStudyModule &module)
Insert given RooMCStudy add-on module to the processing chain of this MCStudy object.
RooArgSet _fitInitParams
List of initial values of fit parameters.
Definition RooMCStudy.h:121
RooPlot * plotParam(const RooRealVar &param, const RooCmdArg &arg1={}, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={})
Plot the distribution of the fitted value of the given parameter on a newly created frame.
std::unique_ptr< RooRealVar > _ngenVar
Definition RooMCStudy.h:124
RooPlot * plotError(const RooRealVar &param, const RooCmdArg &arg1={}, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={})
Plot the distribution of the fit errors for the specified parameter on a newly created frame.
static RooMsgService & instance()
Return reference to singleton instance.
void setGlobalKillBelow(RooFit::MsgLevel level)
RooFit::MsgLevel globalKillBelow() const
Plot frame and a container for graphics objects within that frame.
Definition RooPlot.h:43
static RooPlot * frame(const RooAbsRealLValue &var, double xmin, double xmax, Int_t nBins)
Create a new frame for a given variable in x.
Definition RooPlot.cxx:237
RooAbsRealLValue * getPlotVar() const
Definition RooPlot.h:139
void createInternalPlotVarClone()
Replaces the pointer to the plot variable with a pointer to a clone of the plot variable that is owne...
Definition RooPlot.cxx:1449
Represents the pull of a measurement w.r.t.
Definition RooPullVar.h:24
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
Definition RooRandom.cxx:51
Variable that can be changed from the outside.
Definition RooRealVar.h:37
RooErrorVar * errorVar() const
Return a RooAbsRealLValue representing the error associated with this variable.
TString * format(const RooCmdArg &formatArg) const
Format contents of RooRealVar for pretty printing on RooPlot parameter boxes.
void setBins(Int_t nBins, const char *name=nullptr)
Create a uniform binning under name 'name' for this variable.
Persistable container for RooFit projects.
RooAbsPdf * pdf(RooStringView name) const
Retrieve p.d.f (RooAbsPdf) with given name. A null pointer is returned if not found.
bool import(const RooAbsArg &arg, const RooCmdArg &arg1={}, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}, const RooCmdArg &arg9={})
Import a RooAbsArg object, e.g.
RooFactoryWSTool & factory()
Return instance to factory tool.
RooRealVar * var(RooStringView name) const
Retrieve real-valued variable (RooRealVar) with given name. A null pointer is returned if not found.
virtual Int_t GetSize() const
Return the capacity of the collection, i.e.
A doubly linked list.
Definition TList.h:38
void RecursiveRemove(TObject *obj) override
Remove object from this collection and recursively remove the object from all other objects (and coll...
Definition TList.cxx:764
void Add(TObject *obj) override
Definition TList.h:81
void Delete(Option_t *option="") override
Remove all objects from the list AND delete all heap based objects.
Definition TList.cxx:470
TObject * At(Int_t idx) const override
Returns the object at position idx. Returns 0 if idx is out of range.
Definition TList.cxx:357
The TNamed class is the base class for all named ROOT classes.
Definition TNamed.h:29
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
const char * GetTitle() const override
Returns title of object.
Definition TNamed.h:48
Mother of all ROOT objects.
Definition TObject.h:41
virtual ULong64_t Poisson(Double_t mean)
Generates a random integer N according to a Poisson law.
Definition TRandom.cxx:404
Basic string class.
Definition TString.h:139
TString & Append(const char *cs)
Definition TString.h:574
RooCmdArg AutoRange(const RooAbsData &data, double marginFactor=0.1)
RooCmdArg Label(const char *str)
RooCmdArg AutoSymRange(const RooAbsData &data, double marginFactor=0.1)
RooCmdArg Bins(Int_t nbin)
RooCmdArg Layout(double xmin, double xmax=0.99, double ymin=0.95)
RooCmdArg Constrain(const RooArgSet &params)
RooCmdArg Save(bool flag=true)
RooCmdArg ExternalConstraints(const RooArgSet &constraintPdfs)
RooCmdArg Minos(bool flag=true)
RooCmdArg PrintLevel(Int_t code)
RooCmdArg ConditionalObservables(Args_t &&... argsOrArgSet)
Create a RooCmdArg to declare conditional observables.
RooCmdArg Range(const char *rangeName, bool adjustNorm=true)
MsgLevel
Verbosity level for RooMsgService::StreamConfig in RooMsgService.
OwningPtr< T > makeOwningPtr(std::unique_ptr< T > &&ptr)
Internal helper to turn a std::unique_ptr<T> into an OwningPtr.
Definition Config.h:40
T * OwningPtr
An alias for raw pointers for indicating that the return type of a RooFit function is an owning point...
Definition Config.h:35