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> ProtoData(const RooDataSet&, bool randOrder)
100 <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
101 cycle to protect against systematic biases if the number of generated events does not exactly match the number of events in the prototype dataset
102 at the cost of reduced precision with mu equal to the specified number of events
103</table>
104*/
105RooMCStudy::RooMCStudy(const RooAbsPdf& model, const RooArgSet& observables,
106 const RooCmdArg& arg1, const RooCmdArg& arg2,
107 const RooCmdArg& arg3,const RooCmdArg& arg4,const RooCmdArg& arg5,
108 const RooCmdArg& arg6,const RooCmdArg& arg7,const RooCmdArg& arg8) : TNamed("mcstudy","mcstudy")
109
110{
111 // Stuff all arguments in a list
112 RooLinkedList cmdList;
113 cmdList.Add(const_cast<RooCmdArg*>(&arg1)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg2)) ;
114 cmdList.Add(const_cast<RooCmdArg*>(&arg3)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg4)) ;
115 cmdList.Add(const_cast<RooCmdArg*>(&arg5)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg6)) ;
116 cmdList.Add(const_cast<RooCmdArg*>(&arg7)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg8)) ;
117
118 // Select the pdf-specific commands
119 RooCmdConfig pc("RooMCStudy::RooMCStudy(" + std::string(model.GetName()) + ")");
120
121 pc.defineObject("fitModel","FitModel",0,nullptr) ;
122 pc.defineSet("condObs","ProjectedObservables",0,nullptr) ;
123 pc.defineObject("protoData","PrototypeData",0,nullptr) ;
124 pc.defineSet("cPars","Constrain",0,nullptr) ;
125 pc.defineSet("extCons","ExternalConstraints",0,nullptr) ;
126 pc.defineInt("silence","Silence",0,0) ;
127 pc.defineInt("randProtoData","PrototypeData",0,0) ;
128 pc.defineInt("verboseGen","Verbose",0,0) ;
129 pc.defineInt("extendedGen","Extended",0,0) ;
130 pc.defineInt("binGenData","Binned",0,0) ;
131 pc.defineInt("dummy","FitOptArgs",0,0) ;
132
133 // Process and check varargs
134 pc.process(cmdList) ;
135 if (!pc.ok(true)) {
136 // WVE do something here
137 throw std::string("RooMCStudy::RooMCStudy() Error in parsing arguments passed to constructor") ;
138 return ;
139 }
140
141 // Save fit command options
142 if (pc.hasProcessed("FitOptArgs")) {
143 RooCmdArg* fitOptArg = static_cast<RooCmdArg*>(cmdList.FindObject("FitOptArgs")) ;
144 for (int i=0 ; i<fitOptArg->subArgs().GetSize() ;i++) {
145 _fitOptList.Add(new RooCmdArg(static_cast<RooCmdArg&>(*fitOptArg->subArgs().At(i)))) ;
146 }
147 }
148
149 // Decode command line arguments
150 _silence = pc.getInt("silence") ;
151 _verboseGen = pc.getInt("verboseGen") ;
152 _extendedGen = pc.getInt("extendedGen") ;
153 _binGenData = pc.getInt("binGenData") ;
154 _randProto = pc.getInt("randProtoData") ;
155
156 // Process constraints specifications
157 const RooArgSet* cParsTmp = pc.getSet("cPars") ;
158 const RooArgSet* extCons = pc.getSet("extCons") ;
159
160 auto cPars = std::make_unique<RooArgSet>();
161 if (cParsTmp) {
162 cPars->add(*cParsTmp) ;
163 }
164
165 // If constraints are specified, add to fit options
166 if (cPars) {
168 }
169 if (extCons) {
171 }
172
173 // Make list of all constraints
174 RooArgSet allConstraints ;
175 RooArgSet consPars ;
176 if (cPars) {
177 if (std::unique_ptr<RooArgSet> constraints{model.getAllConstraints(observables,*cPars,true)}) {
178 allConstraints.add(*constraints) ;
179 }
180 }
181
182 // Construct constraint p.d.f
183 if (!allConstraints.empty()) {
184 _constrPdf = std::make_unique<RooProdPdf>("mcs_constr_prod","RooMCStudy constraints product",allConstraints);
185
186 if (cPars) {
187 consPars.add(*cPars) ;
188 } else {
189 RooArgSet params;
190 model.getParameters(&observables, params);
191 RooArgSet cparams;
192 _constrPdf->getObservables(&params, cparams);
193 consPars.add(cparams) ;
194 }
195 _constrGenContext.reset(_constrPdf->genContext(consPars,nullptr,nullptr,_verboseGen));
196
197 _perExptGenParams = true ;
198
199 coutI(Generation) << "RooMCStudy::RooMCStudy: INFO have pdf with constraints, will generate parameters from constraint pdf for each experiment" << std::endl ;
200 }
201
202
203 // Extract generator and fit models
204 _genModel = const_cast<RooAbsPdf*>(&model) ;
205 RooAbsPdf* fitModel = static_cast<RooAbsPdf*>(pc.getObject("fitModel",nullptr)) ;
206 _fitModel = fitModel ? fitModel : _genModel ;
207
208 // Extract conditional observables and prototype data
209 _genProtoData = static_cast<RooDataSet*>(pc.getObject("protoData",nullptr)) ;
210 if (auto condObs = pc.getSet("condObs",nullptr)) {
211 _projDeps.add(*condObs);
212 }
213
214 _dependents.add(observables) ;
215
217 _canAddFitResults = true ;
218
220 oocoutW(_fitModel,Generation) << "RooMCStudy::RooMCStudy: WARNING Using generator option 'e' (Poisson distribution of #events) together " << std::endl
221 << " with a prototype dataset implies incomplete sampling or oversampling of proto data." << std::endl
222 << " Use option \"r\" to randomize prototype dataset order and thus to randomize" << std::endl
223 << " the set of over/undersampled prototype events for each generation cycle." << std::endl ;
224 }
225
227 if (!_binGenData) {
229 _genContext->attach(_genParams) ;
230 }
231
233
234 // Store list of parameters and save initial values separately
237
239
240 // Place holder for NLL
241 _nllVar = std::make_unique<RooRealVar>("NLL","-log(Likelihood)",0);
242
243 // Place holder for number of generated events
244 _ngenVar = std::make_unique<RooRealVar>("ngen","number of generated events",0);
245
246 // Create data set containing parameter values, errors and pulls
247 RooArgSet tmp2(_fitParams) ;
248 tmp2.add(*_nllVar) ;
249 tmp2.add(*_ngenVar) ;
250
251 // Mark all variable to store their errors in the dataset
252 tmp2.setAttribAll("StoreError",true) ;
253 tmp2.setAttribAll("StoreAsymError",true) ;
254 std::string fpdName;
255 if (_fitModel==_genModel) {
256 fpdName = "fitParData_" + std::string(_fitModel->GetName());
257 } else {
258 fpdName= "fitParData_" + std::string(_fitModel->GetName()) + "_" + std::string(_genModel->GetName());
259 }
260
261 _fitParData = std::make_unique<RooDataSet>(fpdName,"Fit Parameters DataSet",tmp2);
262 tmp2.setAttribAll("StoreError",false) ;
263 tmp2.setAttribAll("StoreAsymError",false) ;
264
265 if (_perExptGenParams) {
266 _genParData = std::make_unique<RooDataSet>("genParData","Generated Parameters dataset",_genParams);
267 }
268
269 // Append proto variables to allDependents
270 if (_genProtoData) {
272 }
273
274 // Call module initializers
275 for (auto iter=_modList.begin() ; iter!= _modList.end() ; ++iter) {
276 bool ok = (*iter)->doInitializeInstance(*this) ;
277 if (!ok) {
278 oocoutE(_fitModel,Generation) << "RooMCStudy::ctor: removing study module " << (*iter)->GetName() << " from analysis chain because initialization failed" << std::endl ;
279 iter = _modList.erase(iter) ;
280 }
281 }
282
283}
284
285
286////////////////////////////////////////////////////////////////////////////////
287
289{
292}
293
294
295
296////////////////////////////////////////////////////////////////////////////////
297/// Insert given RooMCStudy add-on module to the processing chain
298/// of this MCStudy object
299
301{
302 module.doInitializeInstance(*this) ;
303 _modList.push_back(&module) ;
304}
305
306
307
308////////////////////////////////////////////////////////////////////////////////
309/// Run engine method. Generate and/or fit, according to flags, 'nSamples' samples of 'nEvtPerSample' events.
310/// If keepGenData is set, all generated data sets will be kept in memory and can be accessed
311/// later via genData().
312///
313/// When generating, data sets will be written out in ascii form if the pattern string is supplied
314/// The pattern, which is a template for snprintf, should look something like "data/toymc_%04d.dat"
315/// and should contain one integer field that encodes the sample serial number.
316///
317/// When fitting only, data sets may optionally be read from ascii files, using the same file
318/// pattern.
319///
320
321bool RooMCStudy::run(bool doGenerate, bool DoFit, Int_t nSamples, Int_t nEvtPerSample, bool keepGenData, const char* asciiFilePat)
322{
324 if (_silence) {
327 }
328
329 for (RooAbsMCStudyModule *mod : _modList) {
330 mod->initializeRun(nSamples) ;
331 }
332
333 int prescale = nSamples>100 ? int(nSamples/100) : 1 ;
334
335 while(nSamples--) {
336
337 if (nSamples%prescale==0) {
338 oocoutP(_fitModel,Generation) << "RooMCStudy::run: " ;
339 if (doGenerate) ooccoutI(_fitModel,Generation) << "Generating " ;
340 if (doGenerate && DoFit) ooccoutI(_fitModel,Generation) << "and " ;
341 if (DoFit) ooccoutI(_fitModel,Generation) << "fitting " ;
342 ooccoutP(_fitModel,Generation) << "sample " << nSamples << std::endl ;
343 }
344
345 std::unique_ptr<RooAbsData> ownedGenSample;
346 _genSample = nullptr;
347 bool existingData = false ;
348 if (doGenerate) {
349 // Generate sample
350 int nEvt(nEvtPerSample) ;
351
352 // Reset generator parameters to initial values
354
355 // If constraints are present, sample generator values from constraints
356 if (_constrPdf) {
357 _genParams.assign(*std::unique_ptr<RooDataSet>{_constrGenContext->generate(1)}->get());
358 }
359
360 // Save generated parameters if required
361 if (_genParData) {
362 _genParData->add(_genParams) ;
363 }
364
365 // Call module before-generation hook
366 for (RooAbsMCStudyModule *mod : _modList) {
367 mod->processBeforeGen(nSamples) ;
368 }
369
370 if (_binGenData) {
371
372 // Calculate the number of (extended) events for this run
373 if (_extendedGen) {
375 nEvt = RooRandom::randomGenerator()->Poisson(nEvtPerSample==0?_nExpGen:nEvtPerSample) ;
376 }
377
378 // Binned generation
379 ownedGenSample = std::unique_ptr<RooDataHist>{_genModel->generateBinned(_dependents,nEvt)};
380
381 } else {
382
383 // Calculate the number of (extended) events for this run
384 if (_extendedGen) {
386 nEvt = RooRandom::randomGenerator()->Poisson(nEvtPerSample==0?_nExpGen:nEvtPerSample) ;
387 }
388
389 // Optional randomization of protodata for this run
391 oocoutI(_fitModel,Generation) << "RooMCStudy: (Re)randomizing event order in prototype dataset (Nevt=" << nEvt << ")" << std::endl ;
393 _genContext->setProtoDataOrder(newOrder) ;
394 delete[] newOrder ;
395 }
396
397 // Actual generation of events
398 if (nEvt>0) {
399 ownedGenSample = std::unique_ptr<RooAbsData>{_genContext->generate(nEvt)};
400 } else {
401 // Make empty dataset
402 ownedGenSample = std::make_unique<RooDataSet>("emptySample","emptySample",_dependents);
403 }
404 }
405
406 _genSample = ownedGenSample.get();
407
408 //} else if (asciiFilePat && &asciiFilePat) { //warning: the address of 'asciiFilePat' will always evaluate as 'true'
409 } else if (asciiFilePat) {
410
411 // Load sample from ASCII file
412 char asciiFile[1024] ;
413 snprintf(asciiFile,1024,asciiFilePat,nSamples) ;
414 RooArgList depList(_allDependents) ;
415 ownedGenSample = std::unique_ptr<RooDataSet>{RooDataSet::read(asciiFile,depList,"q")};
416 _genSample = ownedGenSample.get();
417
418 } else {
419
420 // Load sample from internal list
421 _genSample = static_cast<RooDataSet*>(_genDataList.At(nSamples)) ;
422 existingData = true ;
423 if (!_genSample) {
424 oocoutW(_fitModel,Generation) << "RooMCStudy::run: WARNING: Sample #" << nSamples << " not loaded, skipping" << std::endl ;
425 continue ;
426 }
427 }
428
429 // Save number of generated events
430 _ngenVar->setVal(_genSample->sumEntries()) ;
431
432 // Call module between generation and fitting hook
433 for (RooAbsMCStudyModule *mod : _modList) {
434 mod->processBetweenGenAndFit(nSamples) ;
435 }
436
437 if (DoFit) fitSample(_genSample) ;
438
439 // Call module between generation and fitting hook
440 for (RooAbsMCStudyModule *mod : _modList) {
441 mod->processAfterFit(nSamples) ;
442 }
443
444 // Optionally write to ascii file
445 if (doGenerate && asciiFilePat && *asciiFilePat) {
446 char asciiFile[1024] ;
447 snprintf(asciiFile,1024,asciiFilePat,nSamples) ;
448 if (RooDataSet* unbinnedData = dynamic_cast<RooDataSet*>(_genSample)) {
449 unbinnedData->write(asciiFile) ;
450 } else {
451 coutE(InputArguments) << "RooMCStudy::run(" << GetName() << ") ERROR: ASCII writing of binned datasets is not supported" << std::endl ;
452 }
453 }
454
455 // Add to list or delete
456 if (!existingData) {
457 if (keepGenData) {
458 _genDataList.Add(ownedGenSample.release()) ;
459 }
460 }
461 }
462
463 for (RooAbsMCStudyModule *mod : _modList) {
464 if (RooDataSet* auxData = mod->finalizeRun()) {
465 _fitParData->merge(auxData) ;
466 }
467 }
468
469 _canAddFitResults = false ;
470
471 if (_genParData) {
472 for(RooAbsArg * arg : *_genParData->get()) {
473 _genParData->changeObservableName(arg->GetName(),(std::string(arg->GetName()) + "_gen").c_str());
474 }
475
476 _fitParData->merge(_genParData.get());
477 }
478
479 if (DoFit) calcPulls() ;
480
481 if (_silence) {
483 }
484
485 return false ;
486}
487
488
489
490
491
492
493////////////////////////////////////////////////////////////////////////////////
494/// Generate and fit 'nSamples' samples of 'nEvtPerSample' events.
495/// If keepGenData is set, all generated data sets will be kept in memory and can be accessed
496/// later via genData().
497///
498/// Data sets will be written out in ascii form if the pattern string is supplied.
499/// The pattern, which is a template for snprintf, should look something like "data/toymc_%04d.dat"
500/// and should contain one integer field that encodes the sample serial number.
501///
502
503bool RooMCStudy::generateAndFit(Int_t nSamples, Int_t nEvtPerSample, bool keepGenData, const char* asciiFilePat)
504{
505 // Clear any previous data in memory
506 _fitResList.Delete() ; // even though the fit results are owned by gROOT, we still want to scratch them here.
508 _fitParData->reset() ;
509
510 return run(true,true,nSamples,nEvtPerSample,keepGenData,asciiFilePat) ;
511}
512
513
514
515////////////////////////////////////////////////////////////////////////////////
516/// Generate 'nSamples' samples of 'nEvtPerSample' events.
517/// If keepGenData is set, all generated data sets will be kept in memory
518/// and can be accessed later via genData().
519///
520/// Data sets will be written out in ascii form if the pattern string is supplied.
521/// The pattern, which is a template for snprintf, should look something like "data/toymc_%04d.dat"
522/// and should contain one integer field that encodes the sample serial number.
523///
524
525bool RooMCStudy::generate(Int_t nSamples, Int_t nEvtPerSample, bool keepGenData, const char* asciiFilePat)
526{
527 // Clear any previous data in memory
529
530 return run(true,false,nSamples,nEvtPerSample,keepGenData,asciiFilePat) ;
531}
532
533
534
535////////////////////////////////////////////////////////////////////////////////
536/// Fit 'nSamples' datasets, which are read from ASCII files.
537///
538/// The ascii file pattern, which is a template for snprintf, should look something like "data/toymc_%04d.dat"
539/// and should contain one integer field that encodes the sample serial number.
540///
541
542bool RooMCStudy::fit(Int_t nSamples, const char* asciiFilePat)
543{
544 // Clear any previous data in memory
545 _fitResList.Delete() ; // even though the fit results are owned by gROOT, we still want to scratch them here.
546 _fitParData->reset() ;
547
548 return run(false,true,nSamples,0,false,asciiFilePat) ;
549}
550
551
552
553////////////////////////////////////////////////////////////////////////////////
554/// Fit 'nSamples' datasets, as supplied in 'dataSetList'
555///
556
557bool RooMCStudy::fit(Int_t nSamples, TList& dataSetList)
558{
559 // Clear any previous data in memory
560 _fitResList.Delete() ; // even though the fit results are owned by gROOT, we still want to scratch them here.
562 _fitParData->reset() ;
563
564 // Load list of data sets
565 for(auto * gset : static_range_cast<RooAbsData*>(dataSetList)) {
566 _genDataList.Add(gset) ;
567 }
568
569 return run(false,true,nSamples,0,true,nullptr) ;
570}
571
572
573
574////////////////////////////////////////////////////////////////////////////////
575/// Reset all fit parameters to the initial model
576/// parameters at the time of the RooMCStudy constructor
577
579{
581}
582
583
584
585////////////////////////////////////////////////////////////////////////////////
586/// Internal function. Performs actual fit according to specifications
587
589{
590 // Optionally bin dataset before fitting
591 std::unique_ptr<RooDataHist> ownedDataHist;
593 if (_binGenData) {
594 RooArgSet depList;
595 _fitModel->getObservables(genSample->get(), depList);
596 ownedDataHist = std::make_unique<RooDataHist>(genSample->GetName(),genSample->GetTitle(),depList,*genSample) ;
597 data = ownedDataHist.get();
598 } else {
599 data = genSample ;
600 }
601
602 RooCmdArg save = RooFit::Save() ;
604 RooCmdArg plevel = RooFit::PrintLevel(_silence ? -1 : 1) ;
605
606 RooLinkedList fitOptList(_fitOptList) ;
607 fitOptList.Add(&save) ;
608 if (!_projDeps.empty()) {
609 fitOptList.Add(&condo) ;
610 }
611 fitOptList.Add(&plevel) ;
612 return _fitModel->fitTo(*data,fitOptList);
613}
614
615
616
617////////////////////////////////////////////////////////////////////////////////
618/// Redo fit on 'current' toy sample, or if genSample is not nullptr
619/// do fit on given sample instead
620
622{
623 if (!genSample) {
624 genSample = _genSample ;
625 }
626
627 std::unique_ptr<RooFitResult> fr;
628 if (genSample->sumEntries()>0) {
629 fr = std::unique_ptr<RooFitResult>{doFit(genSample)};
630 }
631
632 return RooFit::makeOwningPtr(std::move(fr));
633}
634
635
636
637////////////////////////////////////////////////////////////////////////////////
638/// Internal method. Fit given dataset with fit model. If fit
639/// converges (TMinuit status code zero) The fit results are appended
640/// to the fit results dataset
641///
642/// If the fit option "r" is supplied, the RooFitResult
643/// objects will always be saved, regardless of the
644/// fit status. RooFitResults objects can be retrieved
645/// later via fitResult().
646///
647
649{
650 // Reset all fit parameters to their initial values
652
653 // Perform actual fit
654 bool ok ;
655 std::unique_ptr<RooFitResult> fr;
656 if (genSample->sumEntries()>0) {
657 fr = std::unique_ptr<RooFitResult>{doFit(genSample)};
658 ok = (fr->status()==0) ;
659 } else {
660 ok = false ;
661 }
662
663 // If fit converged, store parameters and NLL
664 if (ok) {
665 _nllVar->setVal(fr->minNll()) ;
666 RooArgSet tmp(_fitParams) ;
667 tmp.add(*_nllVar) ;
668 tmp.add(*_ngenVar) ;
669
670 _fitParData->add(tmp) ;
671 }
672
673 // Store fit result if requested by user
674 if (_fitOptList.FindObject("Save")) {
675 _fitResList.Add(fr.release()) ;
676 }
677
678 return !ok ;
679}
680
681
682
683////////////////////////////////////////////////////////////////////////////////
684/// Utility function to add fit result from external fit to this RooMCStudy
685/// and process its results through the standard RooMCStudy statistics gathering tools.
686/// This function allows users to run the toy MC generation and/or fitting
687/// in a distributed way and to collect and analyze the results in a RooMCStudy
688/// as if they were run locally.
689///
690/// This method is only functional if this RooMCStudy object is cleanm, i.e. it was not used
691/// to generate and/or fit any samples.
692
694{
695 if (!_canAddFitResults) {
696 oocoutE(_fitModel,InputArguments) << "RooMCStudy::addFitResult: ERROR cannot add fit results in current state" << std::endl ;
697 return true ;
698 }
699
700 // Transfer contents of fit result to fitParams ;
702
703 // If fit converged, store parameters and NLL
704 bool ok = (fr.status()==0) ;
705 if (ok) {
706 _nllVar->setVal(fr.minNll()) ;
707 RooArgSet tmp(_fitParams) ;
708 tmp.add(*_nllVar) ;
709 tmp.add(*_ngenVar) ;
710 _fitParData->add(tmp) ;
711 }
712
713 // Store fit result if requested by user
714 if (_fitOptList.FindObject("Save")) {
715 _fitResList.Add((TObject*)&fr) ;
716 }
717
718 return false ;
719}
720
721
722
723////////////////////////////////////////////////////////////////////////////////
724/// Calculate the pulls for all fit parameters in
725/// the fit results data set, and add them to that dataset.
726
728{
729 for (auto it = _fitParams.begin(); it != _fitParams.end(); ++it) {
730 const auto par = static_cast<RooRealVar*>(*it);
731 _fitParData->addColumn(*std::unique_ptr<RooErrorVar>{par->errorVar()});
732
733 TString name(par->GetName());
734 TString title(par->GetTitle());
735 name.Append("pull") ;
736 title.Append(" Pull") ;
737
738 if (!par->hasError(false)) {
739 coutW(Generation) << "Fit parameter '" << par->GetName() << "' does not have an error."
740 " A pull distribution cannot be generated. This might be caused by the parameter being constant or"
741 " because the fits were not run." << std::endl;
742 continue;
743 }
744
745 // First look in fitParDataset to see if per-experiment generated value has been stored
746 auto genParOrig = static_cast<RooAbsReal*>(_fitParData->get()->find(Form("%s_gen",par->GetName())));
747 if (genParOrig && _perExptGenParams) {
748
749 RooPullVar pull(name,title,*par,*genParOrig) ;
750 _fitParData->addColumn(pull,false) ;
751
752 } else {
753 // If not use fixed generator value
754 genParOrig = static_cast<RooAbsReal*>(_genInitParams.find(par->GetName()));
755
756 if (!genParOrig) {
757 std::size_t index = it - _fitParams.begin();
758 genParOrig = index < _genInitParams.size() ?
759 static_cast<RooAbsReal*>(_genInitParams[index]) :
760 nullptr;
761
762 if (genParOrig) {
763 coutW(Generation) << "The fit parameter '" << par->GetName() << "' is not in the model that was used to generate toy data. "
764 "The parameter '" << genParOrig->GetName() << "'=" << genParOrig->getVal() << " was found at the same position in the generator model."
765 " It will be used to compute pulls."
766 "\nIf this is not desired, the parameters of the generator model need to be renamed or reordered." << std::endl;
767 }
768 }
769
770 if (genParOrig) {
771 std::unique_ptr<RooAbsReal> genPar(static_cast<RooAbsReal*>(genParOrig->Clone("truth")));
772 RooPullVar pull(name,title,*par,*genPar);
773
774 _fitParData->addColumn(pull,false) ;
775 } else {
776 coutE(Generation) << "Cannot generate pull distribution for the fit parameter '" << par->GetName() << "'."
777 "\nNo similar parameter was found in the set of parameters that were used to generate toy data." << std::endl;
778 }
779 }
780 }
781}
782
783
784
785
786////////////////////////////////////////////////////////////////////////////////
787/// Return a RooDataSet containing the post-fit parameters of each toy cycle.
788/// This dataset also contains any additional output that was generated
789/// by study modules that were added to this RooMCStudy.
790/// By default, the two following variables are added (apart from fit parameters):
791/// - NLL: The value of the negative log-likelihood for each run.
792/// - ngen: Number of events generated for each run.
794{
795 if (_canAddFitResults) {
796 calcPulls() ;
797 _canAddFitResults = false ;
798 }
799
800 return *_fitParData ;
801}
802
803
804
805////////////////////////////////////////////////////////////////////////////////
806/// Return an argset with the fit parameters for the given sample number
807///
808/// NB: The fit parameters are only stored for successful fits,
809/// thus the maximum sampleNum can be less that the number
810/// of generated samples and if so, the indices will
811/// be out of synch with genData() and fitResult()
812
813const RooArgSet* RooMCStudy::fitParams(Int_t sampleNum) const
814{
815 // Check if sampleNum is in range
816 if (sampleNum<0 || sampleNum>=_fitParData->numEntries()) {
817 oocoutE(_fitModel,InputArguments) << "RooMCStudy::fitParams: ERROR, invalid sample number: " << sampleNum << std::endl ;
818 return nullptr ;
819 }
820
821 return _fitParData->get(sampleNum) ;
822}
823
824
825
826////////////////////////////////////////////////////////////////////////////////
827/// Return the RooFitResult of the fit with the given run number.
828///
829/// \note Fit results are not saved by default. This requires passing `FitOptions(Save(), ...)`
830/// to the constructor.
832{
833 // Check if sampleNum is in range
834 if (sampleNum<0 || sampleNum>=_fitResList.GetSize()) {
835 oocoutE(_fitModel,InputArguments) << "RooMCStudy::fitResult: ERROR, invalid sample number: " << sampleNum << std::endl ;
836 return nullptr ;
837 }
838
839 // Retrieve fit result object
840 const RooFitResult* fr = static_cast<RooFitResult*>(_fitResList.At(sampleNum)) ;
841 if (fr) {
842 return fr ;
843 } else {
844 oocoutE(_fitModel,InputArguments) << "RooMCStudy::fitResult: ERROR, no fit result saved for sample "
845 << sampleNum << ", did you use the 'r; fit option?" << std::endl ;
846 }
847 return nullptr ;
848}
849
850
851
852////////////////////////////////////////////////////////////////////////////////
853/// Return the given generated dataset. This method will only return datasets
854/// if during the run cycle it was indicated that generator data should be saved.
855
857{
858 // Check that generated data was saved
859 if (_genDataList.GetSize()==0) {
860 oocoutE(_fitModel,InputArguments) << "RooMCStudy::genData() ERROR, generated data was not saved" << std::endl ;
861 return nullptr ;
862 }
863
864 // Check if sampleNum is in range
865 if (sampleNum<0 || sampleNum>=_genDataList.GetSize()) {
866 oocoutE(_fitModel,InputArguments) << "RooMCStudy::genData() ERROR, invalid sample number: " << sampleNum << std::endl ;
867 return nullptr ;
868 }
869
870 return static_cast<RooAbsData*>(_genDataList.At(sampleNum)) ;
871}
872
873
874
875////////////////////////////////////////////////////////////////////////////////
876/// Plot the distribution of fitted values of a parameter. The parameter shown is the one from which the RooPlot
877/// was created, e.g.
878///
879/// RooPlot* frame = param.frame(100,-10,10) ;
880/// mcstudy.paramOn(frame,LineStyle(kDashed)) ;
881///
882/// Any named arguments passed to plotParamOn() are forwarded to the underlying plotOn() call
883
884RooPlot* RooMCStudy::plotParamOn(RooPlot* frame, const RooCmdArg& arg1, const RooCmdArg& arg2, const RooCmdArg& arg3, const RooCmdArg& arg4,
885 const RooCmdArg& arg5, const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8)
886{
887 _fitParData->plotOn(frame,arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8) ;
888 return frame ;
889}
890
891
892
893////////////////////////////////////////////////////////////////////////////////
894/// Plot the distribution of the fitted value of the given parameter on a newly created frame.
895///
896/// <table>
897/// <tr><th> Optional arguments <th>
898/// <tr><td> FrameRange(double lo, double hi) <td> Set range of frame to given specification
899/// <tr><td> FrameBins(int bins) <td> Set default number of bins of frame to given number
900/// <tr><td> Frame() <td> Pass supplied named arguments to RooAbsRealLValue::frame() function. See there
901/// for list of allowed arguments
902/// </table>
903/// If no frame specifications are given, the AutoRange() feature will be used to set the range
904/// Any other named argument is passed to the RooAbsData::plotOn() call. See that function for allowed options
905
906RooPlot* RooMCStudy::plotParam(const char* paramName, const RooCmdArg& arg1, const RooCmdArg& arg2, const RooCmdArg& arg3, const RooCmdArg& arg4,
907 const RooCmdArg& arg5, const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8)
908{
909
910 // Find parameter in fitParDataSet
911 RooRealVar* param = static_cast<RooRealVar*>(_fitParData->get()->find(paramName)) ;
912 if (!param) {
913 oocoutE(_fitModel,InputArguments) << "RooMCStudy::plotParam: ERROR: no parameter defined with name " << paramName << std::endl ;
914 return nullptr ;
915 }
916
917 // Forward to implementation below
918 return plotParam(*param,arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8) ;
919}
920
921
922
923////////////////////////////////////////////////////////////////////////////////
924/// Plot the distribution of the fitted value of the given parameter on a newly created frame.
925/// \copydetails RooMCStudy::plotParam(const char* paramName, const RooCmdArg& arg1, const RooCmdArg& arg2, const RooCmdArg& arg3, const RooCmdArg& arg4,
926/// const RooCmdArg& arg5, const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8)
927
928RooPlot* RooMCStudy::plotParam(const RooRealVar& param, const RooCmdArg& arg1, const RooCmdArg& arg2, const RooCmdArg& arg3, const RooCmdArg& arg4,
929 const RooCmdArg& arg5, const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8)
930{
931 // Stuff all arguments in a list
932 RooLinkedList cmdList;
933 cmdList.Add(const_cast<RooCmdArg*>(&arg1)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg2)) ;
934 cmdList.Add(const_cast<RooCmdArg*>(&arg3)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg4)) ;
935 cmdList.Add(const_cast<RooCmdArg*>(&arg5)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg6)) ;
936 cmdList.Add(const_cast<RooCmdArg*>(&arg7)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg8)) ;
937
938 RooPlot* frame = makeFrameAndPlotCmd(param, cmdList) ;
939 if (frame) {
940 _fitParData->plotOn(frame, cmdList) ;
941 }
942
943 return frame ;
944}
945
946
947
948////////////////////////////////////////////////////////////////////////////////
949/// Plot the distribution of the -log(L) values on a newly created frame.
950///
951/// <table>
952/// <tr><th> Optional arguments <th>
953/// <tr><td> FrameRange(double lo, double hi) <td> Set range of frame to given specification
954/// <tr><td> FrameBins(int bins) <td> Set default number of bins of frame to given number
955/// <tr><td> Frame() <td> Pass supplied named arguments to RooAbsRealLValue::frame() function. See there
956/// for list of allowed arguments
957/// </table>
958///
959/// If no frame specifications are given, the AutoRange() feature will be used to set the range.
960/// Any other named argument is passed to the RooAbsData::plotOn() call. See that function for allowed options
961
963 const RooCmdArg& arg3, const RooCmdArg& arg4,
964 const RooCmdArg& arg5, const RooCmdArg& arg6,
965 const RooCmdArg& arg7, const RooCmdArg& arg8)
966{
967 return plotParam(*_nllVar,arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8) ;
968}
969
970
971
972////////////////////////////////////////////////////////////////////////////////
973/// Plot the distribution of the fit errors for the specified parameter on a newly created frame.
974///
975/// <table>
976/// <tr><th> Optional arguments <th>
977/// <tr><td> FrameRange(double lo, double hi) <td> Set range of frame to given specification
978/// <tr><td> FrameBins(int bins) <td> Set default number of bins of frame to given number
979/// <tr><td> Frame() <td> Pass supplied named arguments to RooAbsRealLValue::frame() function. See there
980/// for list of allowed arguments
981/// </table>
982///
983/// If no frame specifications are given, the AutoRange() feature will be used to set a default range.
984/// Any other named argument is passed to the RooAbsData::plotOn() call. See that function for allowed options.
985
986RooPlot* RooMCStudy::plotError(const RooRealVar& param, const RooCmdArg& arg1, const RooCmdArg& arg2,
987 const RooCmdArg& arg3, const RooCmdArg& arg4,
988 const RooCmdArg& arg5, const RooCmdArg& arg6,
989 const RooCmdArg& arg7, const RooCmdArg& arg8)
990{
991 if (_canAddFitResults) {
992 calcPulls() ;
993 _canAddFitResults=false ;
994 }
995
996 std::unique_ptr<RooErrorVar> evar{param.errorVar()};
997 std::unique_ptr<RooAbsArg> evar_rrv{evar->createFundamental()};
998 RooPlot* frame = plotParam(static_cast<RooRealVar&>(*evar_rrv),arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8) ;
999
1000 // To make sure the frame has no dangling pointer to evar_rrv.
1002
1003 return frame ;
1004}
1005
1006namespace {
1007
1008// Fits a Gaussian to the pull distribution, plots the fit and prints the fit
1009// parameters on the canvas. Implementation detail of RooMCStudy::plotPull().
1010void fitGaussToPulls(RooPlot& frame, RooDataSet& fitParData)
1011{
1012 // Build the Gaussian fit mode for the pulls, then fit it and plot it. We
1013 // have to use the RooWorkspace factory here, because different from the
1014 // RooMCStudy class, the RooGaussian is not in RooFitCore.
1015 RooWorkspace ws;
1016 auto plotVar = frame.getPlotVar();
1017 const std::string plotVarName = plotVar->GetName();
1018 ws.import(*plotVar);
1019 ws.factory("Gaussian::pullGauss(" + plotVarName + ", pullMean[0.0, -10.0, 10.0], pullSigma[1.0, 0.1, 5.0])");
1020
1021 RooRealVar& pullMean = *ws.var("pullMean");
1022 RooRealVar& pullSigma = *ws.var("pullSigma");
1023 RooAbsPdf& pullGauss = *ws.pdf("pullGauss");
1024
1025 pullGauss.fitTo(fitParData, RooFit::Minos(false), RooFit::PrintLevel(-1)) ;
1026 pullGauss.plotOn(&frame) ;
1027
1028 // Instead of using paramOn() without command arguments to plot the fit
1029 // parameters, we are building the parameter label ourselves for more
1030 // flexibility and pass this together with an appropriate layout
1031 // parametrization to paramOn().
1032 const int sigDigits = 2;
1033 const char * options = "ELU";
1034 std::stringstream ss;
1035 ss << "Fit parameters:\n"
1036 << "#mu: " << *std::unique_ptr<TString>{pullMean.format(sigDigits, options)}
1037 << "\n#sigma: " << *std::unique_ptr<TString>{pullSigma.format(sigDigits, options)};
1038 // We set the parameters constant to disable the default label. Still, we
1039 // use param() on as a wrapper for the text box generation.
1040 pullMean.setConstant(true);
1041 pullSigma.setConstant(true);
1042 pullGauss.paramOn(&frame, RooFit::Label(ss.str().c_str()), RooFit::Layout(0.60, 0.9, 0.9));
1043}
1044
1045} // namespace
1046
1047
1048////////////////////////////////////////////////////////////////////////////////
1049/// Plot the distribution of pull values for the specified parameter on a newly created frame. If asymmetric
1050/// errors are calculated in the fit (by MINOS) those will be used in the pull calculation.
1051///
1052/// If the parameters of the models for generation and fit differ, simple heuristics are used to find the
1053/// corresponding parameters:
1054/// - Parameters have the same name: They will be used to compute pulls.
1055/// - Parameters have different names: The position of the fit parameter in the set of fit parameters will be
1056/// computed. The parameter at the same position in the set of generator parameters will be used.
1057///
1058/// Further options:
1059/// <table>
1060/// <tr><th> Arguments <th> Effect
1061/// <tr><td> FrameRange(double lo, double hi) <td> Set range of frame to given specification
1062/// <tr><td> FrameBins(int bins) <td> Set default number of bins of frame to given number
1063/// <tr><td> Frame() <td> Pass supplied named arguments to RooAbsRealLValue::frame() function. See there
1064/// for list of allowed arguments
1065/// <tr><td> FitGauss(bool flag) <td> Add a gaussian fit to the frame
1066/// </table>
1067///
1068/// If no frame specifications are given, the AutoSymRange() feature will be used to set a default range.
1069/// Any other named argument is passed to the RooAbsData::plotOn(). See that function for allowed options.
1070///
1071/// If you want to have more control over the Gaussian fit to the pull
1072/// distribution, you can also do it after the call to plotPull():
1073///
1074/// ~~~ {.cpp}
1075/// RooPlot *frame = mcstudy->plotPull(myVariable, RooFit::Bins(40), RooFit::FitGauss(false));
1076/// RooRealVar pullMean("pullMean","Mean of pull",0,-10,10) ;
1077/// RooRealVar pullSigma("pullSigma","Width of pull",1,0.1,5) ;
1078/// pullMean.setPlotLabel("pull #mu"); // optional (to get nicer plot labels if you want)
1079/// pullSigma.setPlotLabel("pull #sigma"); // optional
1080/// RooGaussian pullGauss("pullGauss","Gaussian of pull", *frame->getPlotVar(), pullMean, pullSigma);
1081/// pullGauss.fitTo(const_cast<RooDataSet&>(mcstudy->fitParDataSet()),
1082/// RooFit::Minos(0), RooFit::PrintLevel(-1)) ;
1083/// pullGauss.plotOn(frame) ;
1084/// pullGauss.paramOn(frame, RooFit::Layout(0.65, 0.9, 0.9)); // optionally specify label position (xmin, xmax, ymax)
1085/// ~~~
1086
1087RooPlot* RooMCStudy::plotPull(const RooRealVar& param, const RooCmdArg& arg1, const RooCmdArg& arg2,
1088 const RooCmdArg& arg3, const RooCmdArg& arg4,
1089 const RooCmdArg& arg5, const RooCmdArg& arg6,
1090 const RooCmdArg& arg7, const RooCmdArg& arg8)
1091{
1092 // Stuff all arguments in a list
1093 RooLinkedList cmdList;
1094 cmdList.Add(const_cast<RooCmdArg*>(&arg1)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg2)) ;
1095 cmdList.Add(const_cast<RooCmdArg*>(&arg3)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg4)) ;
1096 cmdList.Add(const_cast<RooCmdArg*>(&arg5)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg6)) ;
1097 cmdList.Add(const_cast<RooCmdArg*>(&arg7)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg8)) ;
1098
1099 TString name(param.GetName());
1100 TString title(param.GetTitle());
1101 name.Append("pull") ; title.Append(" Pull") ;
1102 RooRealVar pvar(name,title,-100,100) ;
1103 pvar.setBins(100) ;
1104
1105
1106 RooPlot* frame = makeFrameAndPlotCmd(pvar, cmdList, true) ;
1107 if (frame) {
1108
1109 // Pick up optional FitGauss command from list
1110 RooCmdConfig pc("RooMCStudy::plotPull(" + std::string(_genModel->GetName()) + ")");
1111 pc.defineInt("fitGauss","FitGauss",0,0) ;
1112 pc.allowUndefined() ;
1113 pc.process(cmdList) ;
1114 bool fitGauss=pc.getInt("fitGauss") ;
1115
1116 // Pass stripped command list to plotOn()
1117 RooCmdConfig::stripCmdList(cmdList,"FitGauss") ;
1118 const bool success = _fitParData->plotOn(frame,cmdList) ;
1119
1120 if (!success) {
1121 coutF(Plotting) << "No pull distribution for the parameter '" << param.GetName() << "'. Check logs for errors." << std::endl;
1122 return frame;
1123 }
1124
1125 // Add Gaussian fit if requested
1126 if (fitGauss) {
1127 fitGaussToPulls(*frame, *_fitParData);
1128 }
1129
1130 // To make sure the frame has no dangling pointer to pvar.
1132 }
1133 return frame;
1134}
1135
1136
1137
1138////////////////////////////////////////////////////////////////////////////////
1139/// Internal function. Construct RooPlot from given parameter and modify the list of named
1140/// arguments 'cmdList' to only contain the plot arguments that should be forwarded to
1141/// RooAbsData::plotOn()
1142
1143RooPlot* RooMCStudy::makeFrameAndPlotCmd(const RooRealVar& param, RooLinkedList& cmdList, bool symRange) const
1144{
1145 // Select the frame-specific commands
1146 RooCmdConfig pc("RooMCStudy::plotParam(" + std::string(_genModel->GetName()) + ")");
1147 pc.defineInt("nbins","Bins",0,0) ;
1148 pc.defineDouble("xlo","Range",0,0) ;
1149 pc.defineDouble("xhi","Range",1,0) ;
1150 pc.defineInt("dummy","FrameArgs",0,0) ;
1151 pc.defineMutex("Bins","FrameArgs") ;
1152 pc.defineMutex("Range","FrameArgs") ;
1153
1154 // Process and check varargs
1155 pc.allowUndefined() ;
1156 pc.process(cmdList) ;
1157 if (!pc.ok(true)) {
1158 return nullptr ;
1159 }
1160
1161 // Make frame according to specs
1162 Int_t nbins = pc.getInt("nbins") ;
1163 double xlo = pc.getDouble("xlo") ;
1164 double xhi = pc.getDouble("xhi") ;
1165 RooPlot* frame ;
1166
1167 if (pc.hasProcessed("FrameArgs")) {
1168 // Explicit frame arguments are given, pass them on
1169 RooCmdArg* frameArg = static_cast<RooCmdArg*>(cmdList.FindObject("FrameArgs")) ;
1170 frame = param.frame(frameArg->subArgs()) ;
1171 } else {
1172 // FrameBins, FrameRange or none are given, build custom frame command list
1173 RooCmdArg bins = RooFit::Bins(nbins) ;
1174 RooCmdArg range = RooFit::Range(xlo,xhi) ;
1175 RooCmdArg autoRange = symRange ? RooFit::AutoSymRange(*_fitParData,0.2) : RooFit::AutoRange(*_fitParData,0.2) ;
1176 RooLinkedList frameCmdList ;
1177
1178 if (pc.hasProcessed("Bins")) frameCmdList.Add(&bins) ;
1179 if (pc.hasProcessed("Range")) {
1180 frameCmdList.Add(&range) ;
1181 } else {
1182 frameCmdList.Add(&autoRange) ;
1183 }
1184 frame = param.frame(frameCmdList) ;
1185 }
1186
1187 // Filter frame command from list and pass on to plotOn()
1188 RooCmdConfig::stripCmdList(cmdList,"FrameArgs,Bins,Range") ;
1189
1190 return frame ;
1191}
1192
1193
1194
1195////////////////////////////////////////////////////////////////////////////////
1196/// Create a RooPlot of the -log(L) distribution in the range lo-hi
1197/// with 'nBins' bins
1198
1199RooPlot* RooMCStudy::plotNLL(double lo, double hi, Int_t nBins)
1200{
1201 RooPlot* frame = _nllVar->frame(lo,hi,nBins) ;
1202
1203 _fitParData->plotOn(frame) ;
1204 return frame ;
1205}
1206
1207
1208
1209////////////////////////////////////////////////////////////////////////////////
1210/// Create a RooPlot of the distribution of the fitted errors of the given parameter.
1211/// The frame is created with a range [lo,hi] and plotted data will be binned in 'nbins' bins
1212
1213RooPlot* RooMCStudy::plotError(const RooRealVar& param, double lo, double hi, Int_t nbins)
1214{
1215 if (_canAddFitResults) {
1216 calcPulls() ;
1217 _canAddFitResults=false ;
1218 }
1219
1220 std::unique_ptr<RooErrorVar> evar{param.errorVar()};
1221 RooPlot* frame = evar->frame(lo,hi,nbins) ;
1222 _fitParData->plotOn(frame) ;
1223
1224 return frame ;
1225}
1226
1227
1228
1229////////////////////////////////////////////////////////////////////////////////
1230/// Create a RooPlot of the pull distribution for the given
1231/// parameter. The range lo-hi is plotted in nbins. If fitGauss is
1232/// set, an unbinned ML fit of the distribution to a Gaussian p.d.f
1233/// is performed. The fit result is overlaid on the returned RooPlot
1234/// and a box with the fitted mean and sigma is added.
1235///
1236/// If the parameters of the models for generation and fit differ, simple heuristics are used to find the
1237/// corresponding parameters:
1238/// - Parameters have the same name: They will be used to compute pulls.
1239/// - Parameters have different names: The position of the fit parameter in the set of fit parameters will be
1240/// computed. The parameter at the same position in the set of generator parameters will be used.
1241
1242RooPlot* RooMCStudy::plotPull(const RooRealVar& param, double lo, double hi, Int_t nbins, bool fitGauss)
1243{
1244 if (_canAddFitResults) {
1245 calcPulls() ;
1246 _canAddFitResults=false ;
1247 }
1248
1249 TString name(param.GetName());
1250 TString title(param.GetTitle());
1251 name.Append("pull") ; title.Append(" Pull") ;
1252 RooRealVar pvar(name,title,lo,hi) ;
1253 pvar.setBins(nbins) ;
1254
1255 RooPlot* frame = pvar.frame() ;
1256 const bool success = _fitParData->plotOn(frame);
1257
1258 if (!success) {
1259 coutF(Plotting) << "No pull distribution for the parameter '" << param.GetName() << "'. Check logs for errors." << std::endl;
1260 return frame;
1261 }
1262
1263 if (fitGauss) {
1264 fitGaussToPulls(*frame, *_fitParData);
1265 }
1266
1267 return frame ;
1268}
1269
1270
1271////////////////////////////////////////////////////////////////////////////////
1272/// If one of the TObject we have a referenced to is deleted, remove the
1273/// reference.
1274
1276{
1280 if (_ngenVar.get() == obj) _ngenVar.reset();
1281
1282 if (_fitParData) _fitParData->RecursiveRemove(obj);
1283 if (_fitParData.get() == obj) _fitParData.reset();
1284
1285 if (_genParData) _genParData->RecursiveRemove(obj);
1286 if (_genParData.get() == obj) _genParData.reset();
1287}
1288
#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:48
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:762
void Add(TObject *obj) override
Definition TList.h:83
void Delete(Option_t *option="") override
Remove all objects from the list AND delete all heap based objects.
Definition TList.cxx:468
TObject * At(Int_t idx) const override
Returns the object at position idx. Returns 0 if idx is out of range.
Definition TList.cxx:355
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:572
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