Logo ROOT  
Reference Guide
RooBarlowBeestonLL.cxx
Go to the documentation of this file.
1// @(#)root/roostats:$Id$
2// Author: Kyle Cranmer, George Lewis
3/*************************************************************************
4 * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
5 * All rights reserved. *
6 * *
7 * For the licensing terms see $ROOTSYS/LICENSE. *
8 * For the list of contributors see $ROOTSYS/README/CREDITS. *
9 *************************************************************************/
10
11//////////////////////////////////////////////////////////////////////////////
12/** \class RooStats::HistFactory::RooBarlowBeestonLL
13 * \ingroup HistFactory
14//
15// Class RooBarlowBeestonLL implements the profile likelihood estimator for
16// a given likelihood and set of parameters of interest. The value return by
17// RooBarlowBeestonLL is the input likelihood nll minimized w.r.t all nuisance parameters
18// (which are all parameters except for those listed in the constructor) minus
19// the -log(L) of the best fit. Note that this function is slow to evaluate
20// as a MIGRAD minimization step is executed for each function evaluation
21*/
22
23#include <stdexcept>
24#include <cmath>
25#include <iostream>
26
28#include "RooAbsReal.h"
29#include "RooMsgService.h"
30#include "RooRealVar.h"
31
33#include "RooCategory.h"
34#include "RooSimultaneous.h"
35#include "RooArgList.h"
36
39
40#include <TMath.h>
41
42using namespace std ;
43
45
46
47////////////////////////////////////////////////////////////////////////////////
48
50 RooAbsReal("RooBarlowBeestonLL","RooBarlowBeestonLL"),
51 _nll(),
52// _obs("paramOfInterest","Parameters of interest",this),
53// _par("nuisanceParam","Nuisance parameters",this,false,false),
54 _pdf(nullptr), _data(nullptr)
55{
56 // Default constructor
57 // Should only be used by proof.
58}
59
60
61////////////////////////////////////////////////////////////////////////////////
62
64 RooAbsReal& nllIn /*, const RooArgSet& observables*/) :
65 RooAbsReal(name,title),
66 _nll("input","-log(L) function",this,nllIn),
67 // _obs("paramOfInterest","Parameters of interest",this),
68 // _par("nuisanceParam","Nuisance parameters",this,false,false),
69 _pdf(nullptr), _data(nullptr)
70{
71 // Constructor of profile likelihood given input likelihood nll w.r.t
72 // the given set of variables. The input log likelihood is minimized w.r.t
73 // to all other variables of the likelihood at each evaluation and the
74 // value of the global log likelihood minimum is always subtracted.
75
76 // Determine actual parameters and observables
77 /*
78 RooArgSet* actualObs = nllIn.getObservables(observables) ;
79 RooArgSet* actualPars = nllIn.getParameters(observables) ;
80
81 _obs.add(*actualObs) ;
82 _par.add(*actualPars) ;
83
84 delete actualObs ;
85 delete actualPars ;
86 */
87}
88
89
90
91////////////////////////////////////////////////////////////////////////////////
92
94 RooAbsReal(other,name),
95 _nll("nll",this,other._nll),
96 // _obs("obs",this,other._obs),
97 // _par("par",this,other._par),
98 _pdf(nullptr), _data(nullptr),
99 _paramFixed(other._paramFixed)
100{
101 // Copy constructor
102
103 // _paramAbsMin.addClone(other._paramAbsMin) ;
104 // _obsAbsMin.addClone(other._obsAbsMin) ;
105
106}
107
108
109
110////////////////////////////////////////////////////////////////////////////////
111/// Destructor
112
114{
115 // Delete instance of minuit if it was ever instantiated
116 // if (_minuit) {
117 // delete _minuit ;
118 // }
119
120
121 //delete _piter ;
122 //delete _oiter ;
123}
124
125
126////////////////////////////////////////////////////////////////////////////////
127
129 for (auto const *var : static_range_cast<RooRealVar *>(*bin_center)) {
130 RooRealVar* target = (RooRealVar*) observables->find(var->GetName()) ;
131 target->setVal(var->getVal()) ;
132 }
133}
134
135
136////////////////////////////////////////////////////////////////////////////////
137
139 bool verbose=false;
140
141 if(!_data) {
142 std::cout << "Error: Must initialize data before initializing cache" << std::endl;
143 throw std::runtime_error("Uninitialized Data");
144 }
145 if(!_pdf) {
146 std::cout << "Error: Must initialize model pdf before initializing cache" << std::endl;
147 throw std::runtime_error("Uninitialized model pdf");
148 }
149
150 // Get the data bin values for all channels and bins
151 std::map< std::string, std::vector<double> > ChannelBinDataMap;
152 getDataValuesForObservables( ChannelBinDataMap, _data, _pdf );
153
154 // Get a list of constraint terms
155 RooArgList obsTerms;
156 RooArgList constraints;
157 RooArgSet* obsSet = _pdf->getObservables(*_data);
158 FactorizeHistFactoryPdf(*obsSet, *_pdf, obsTerms, constraints);
159
160 if( obsTerms.empty() ) {
161 std::cout << "Error: Found no observable terms with pdf: " << _pdf->GetName()
162 << " using dataset: " << _data->GetName() << std::endl;
163 return;
164 }
165 if( constraints.empty() ) {
166 std::cout << "Error: Found no constraint terms with pdf: " << _pdf->GetName()
167 << " using dataset: " << _data->GetName() << std::endl;
168 return;
169 }
170
171 /*
172 // Get the channels for this pdf
173 RooArgSet* channels = new RooArgSet();
174 RooArgSet* channelsWithConstraints = new RooArgSet();
175 getChannelsFromModel( _pdf, channels, channelsWithConstraints );
176 */
177
178 // Loop over the channels
179 RooSimultaneous* simPdf = (RooSimultaneous*) _pdf;
180 RooCategory* channelCat = (RooCategory*) (&simPdf->indexCat());
181 for (const auto& nameIdx : *channelCat) {
182
183 // Warning: channel cat name is not necesarily the same name
184 // as the pdf's (for example, if someone does edits)
185 RooAbsPdf* channelPdf = simPdf->getPdf(nameIdx.first.c_str());
186 std::string channel_name = channelPdf->GetName();
187
188 // First, we check if this channel uses Stat Uncertainties:
189 RooArgList* gammas = new RooArgList();
190 ParamHistFunc* param_func=nullptr;
191 bool hasStatUncert = getStatUncertaintyFromChannel( channelPdf, param_func, gammas );
192 if( ! hasStatUncert ) {
193 if(verbose) {
194 std::cout << "Channel: " << channel_name
195 << " doesn't have statistical uncertainties"
196 << std::endl;
197 }
198 continue;
199 }
200 else {
201 if(verbose) std::cout << "Found ParamHistFunc: " << param_func->GetName() << std::endl;
202 }
203
204 // Now, loop over the bins in this channel
205 // To Do: Check that the index convention
206 // still works for 2-d (ie matches the
207 // convention in ParamHistFunc, etc)
208 int num_bins = param_func->numBins();
209
210 // Initialize the vector to the number of bins, allowing
211 // us to skip gamma's if they are constant
212
213 std::vector<BarlowCache> temp_cache( num_bins );
214 bool channel_has_stat_uncertainty=false;
215
216 for( Int_t bin_index = 0; bin_index < num_bins; ++bin_index ) {
217
218 // Create a cache object
219 BarlowCache cache;
220
221 // Get the gamma for this bin, and skip if it's constant
222 RooRealVar* gamma_stat = dynamic_cast<RooRealVar*>(&(param_func->getParameter(bin_index)));
223 if(!gamma_stat) throw std::runtime_error("ParamHistFunc contains non-RooRealVar, not supported in RooBarlowBeestonLL");
224 if( gamma_stat->isConstant() ) {
225 if(verbose) std::cout << "Ignoring constant gamma: " << gamma_stat->GetName() << std::endl;
226 continue;
227 }
228 else {
229 cache.hasStatUncert=true;
230 channel_has_stat_uncertainty=true;
231 cache.gamma = gamma_stat;
232 _statUncertParams.insert( gamma_stat->GetName() );
233 }
234
235 // Store a snapshot of the bin center
236 RooArgSet* bin_center = (RooArgSet*) param_func->get( bin_index )->snapshot();
237 cache.bin_center = bin_center;
238 cache.observables = obsSet;
239
240 cache.binVolume = param_func->binVolume();
241
242 // Delete this part, simply a comment
243 RooArgList obs_list( *(cache.bin_center) );
244
245 // Get the gamma's constraint term
246 RooAbsReal* pois_mean = nullptr;
247 RooRealVar* tau = nullptr;
248 getStatUncertaintyConstraintTerm( &constraints, gamma_stat, pois_mean, tau );
249 if( !tau || !pois_mean ) {
250 std::cout << "Failed to find pois mean or tau parameter for " << gamma_stat->GetName() << std::endl;
251 }
252 else {
253 if(verbose) std::cout << "Found pois mean and tau for parameter: " << gamma_stat->GetName()
254 << " tau: " << tau->GetName() << " " << tau->getVal()
255 << " pois_mean: " << pois_mean->GetName() << " " << pois_mean->getVal()
256 << std::endl;
257 }
258
259 cache.tau = tau;
260 cache.nom_pois_mean = pois_mean;
261
262 // Get the RooRealSumPdf
263 RooAbsPdf* sum_pdf = getSumPdfFromChannel( channelPdf );
264 if( sum_pdf == nullptr ) {
265 std::cout << "Failed to find RooRealSumPdf in channel " << channel_name
266 << ", therefor skipping this channel for analytic uncertainty minimization"
267 << std::endl;
268 channel_has_stat_uncertainty=false;
269 break;
270 }
271 cache.sumPdf = sum_pdf;
272
273 // And set the data value for this bin
274 if( ChannelBinDataMap.find(channel_name) == ChannelBinDataMap.end() ) {
275 std::cout << "Error: channel with name: " << channel_name
276 << " not found in BinDataMap" << std::endl;
277 throw runtime_error("BinDataMap");
278 }
279 double nData = ChannelBinDataMap[channel_name].at(bin_index);
280 cache.nData = nData;
281
282 temp_cache.at(bin_index) = cache;
283 //_barlowCache[channel_name].at(bin_index) = cache;
284
285 } // End: Loop over bins
286
287 if( channel_has_stat_uncertainty ) {
288 std::cout << "Adding channel: " << channel_name
289 << " to the barlow cache" << std::endl;
290 _barlowCache[channel_name] = temp_cache;
291 }
292
293
294 } // End: Loop over channels
295
296
297
298 // Successfully initialized the cache
299 // Printing some info
300 /*
301 std::map< std::string, std::vector< BarlowCache > >::iterator iter_cache;
302 for( iter_cache = _barlowCache.begin(); iter_cache != _barlowCache.end(); ++iter_cache ) {
303
304 std::string channel_name = (*iter_cache).first;
305 std::vector< BarlowCache >& channel_cache = (*iter_cache).second;
306
307
308 for( unsigned int i = 0; i < channel_cache.size(); ++i ) {
309 BarlowCache& bin_cache = channel_cache.at(i);
310
311
312 RooRealVar* gamma = bin_cache.gamma;
313 RooRealVar* tau = bin_cache.tau;
314 RooAbsReal* pois_mean = bin_cache.nom_pois_mean;
315 RooAbsPdf* sum_pdf = (RooAbsPdf*) bin_cache.sumPdf;
316 double binVolume = bin_cache.binVolume;
317
318
319 if( !bin_cache.hasStatUncert ) {
320 std::cout << "Barlow Cache for Channel: " << channel_name
321 << " Bin: " << i
322 << " NO STAT UNCERT"
323 << std::endl;
324 }
325 else {
326 std::cout << "Barlow Cache for Channel: " << channel_name
327 << " Bin: " << i
328 << " gamma: " << gamma->GetName()
329 << " tau: " << tau->GetName()
330 << " pois_mean: " << pois_mean->GetName()
331 << " sum_pdf: " << sum_pdf->GetName()
332 << " binVolume: " << binVolume
333 << std::endl;
334 }
335
336 }
337 }
338 */
339
340}
341
342
343////////////////////////////////////////////////////////////////////////////////
344
346 RooArgSet& outputSet,
347 bool stripDisconnected) const {
348 bool errorInBaseCall = RooAbsArg::getParameters( depList, outputSet, stripDisconnected );
349
350 RooArgSet toRemove;
351 toRemove.reserve( _statUncertParams.size());
352
353 for (auto const& arg : outputSet) {
354
355 // If there is a gamma in the name,
356 // strip it from the list of dependencies
357
358 if( _statUncertParams.find(arg->GetName()) != _statUncertParams.end() ) {
359 toRemove.add( *arg );
360 }
361 }
362
363 for( auto& arg : toRemove) outputSet.remove( *arg, true );
364
365 return errorInBaseCall || false;
366
367}
368
369
370/*
371////////////////////////////////////////////////////////////////////////////////
372
373const RooArgSet& RooStats::HistFactory::RooBarlowBeestonLL::bestFitParams() const
374{
375 validateAbsMin() ;
376 return _paramAbsMin ;
377}
378
379
380////////////////////////////////////////////////////////////////////////////////
381
382const RooArgSet& RooStats::HistFactory::RooBarlowBeestonLL::bestFitObs() const
383{
384 validateAbsMin() ;
385 return _obsAbsMin ;
386}
387*/
388
389
390
391////////////////////////////////////////////////////////////////////////////////
392/// Optimized implementation of createProfile for profile likelihoods.
393/// Return profile of original function in terms of stated parameters
394/// of interest rather than profiling recursively.
395
396/*
397RooAbsReal* RooStats::HistFactory::RooBarlowBeestonLL::createProfile(const RooArgSet& paramsOfInterest)
398{
399 return nll().createProfile(paramsOfInterest) ;
400}
401*/
402
403
404/*
405void RooStats::HistFactory::RooBarlowBeestonLL::FactorizePdf(const RooArgSet &observables, RooAbsPdf &pdf, RooArgList &obsTerms, RooArgList &constraints) const {
406 // utility function to factorize constraint terms from a pdf
407 // (from G. Petrucciani)
408 const std::type_info & id = typeid(pdf);
409 if (id == typeid(RooProdPdf)) {
410 RooProdPdf *prod = dynamic_cast<RooProdPdf *>(&pdf);
411 RooArgList list(prod->pdfList());
412 for (int i = 0, n = list.getSize(); i < n; ++i) {
413 RooAbsPdf *pdfi = (RooAbsPdf *) list.at(i);
414 FactorizePdf(observables, *pdfi, obsTerms, constraints);
415 }
416 } else if (id == typeid(RooSimultaneous) ) { //|| id == typeid(RooSimultaneousOpt)) {
417 RooSimultaneous *sim = dynamic_cast<RooSimultaneous *>(&pdf);
418 RooAbsCategoryLValue *cat = (RooAbsCategoryLValue *) sim->indexCat().Clone();
419 for (int ic = 0, nc = cat->numBins((const char *)0); ic < nc; ++ic) {
420 cat->setBin(ic);
421 FactorizePdf(observables, *sim->getPdf(cat->getLabel()), obsTerms, constraints);
422 }
423 delete cat;
424 } else if (pdf.dependsOn(observables)) {
425 if (!obsTerms.contains(pdf)) obsTerms.add(pdf);
426 } else {
427 if (!constraints.contains(pdf)) constraints.add(pdf);
428 }
429}
430*/
431
432
433
434////////////////////////////////////////////////////////////////////////////////
435
437{
438 // Loop over the channels (keys to the map)
439 //clock_t time_before_setVal, time_after_setVal;
440 //time_before_setVal=clock();
441 std::map< std::string, std::vector< BarlowCache > >::iterator iter_cache;
442 for( iter_cache = _barlowCache.begin(); iter_cache != _barlowCache.end(); ++iter_cache ) {
443
444 std::string channel_name = (*iter_cache).first;
445 std::vector< BarlowCache >& channel_cache = (*iter_cache).second;
446
447 /* Slower way to find the channel vector:
448 // Get the vector of bin uncertainty caches for this channel
449 if( _barlowCache.find( channel_name ) == _barlowCache.end() ) {
450 std::cout << "Error: channel: " << channel_name
451 << " not found in barlow Cache" << std::endl;
452 throw runtime_error("Channel not in barlow cache");
453 }
454
455 std::vector< BarlowCache >& channel_cache = _barlowCache[ channel_name ];
456 */
457
458 // Loop over the bins in the cache
459 // Set all gamma's to 0
460 for( unsigned int i = 0; i < channel_cache.size(); ++i ) {
461 BarlowCache& bin_cache = channel_cache.at(i);
462 if( !bin_cache.hasStatUncert ) continue;
463 RooRealVar* gamma = bin_cache.gamma;
464 gamma->setVal(0.0);
465 }
466 std::vector< double > nu_b_vec( channel_cache.size() );
467 for( unsigned int i = 0; i < channel_cache.size(); ++i ) {
468 BarlowCache& bin_cache = channel_cache.at(i);
469 if( !bin_cache.hasStatUncert ) continue;
470
471 RooAbsPdf* sum_pdf = (RooAbsPdf*) bin_cache.sumPdf;
472 RooArgSet* obsSet = bin_cache.observables;
473 double binVolume = bin_cache.binVolume;
474
475 bin_cache.SetBinCenter();
476 double nu_b = sum_pdf->getVal(*obsSet)*sum_pdf->expectedEvents(*obsSet)*binVolume;
477 nu_b_vec.at(i) = nu_b;
478 }
479
480 // Loop over the bins in the cache
481 // Set all gamma's to 1
482 for( unsigned int i = 0; i < channel_cache.size(); ++i ) {
483 BarlowCache& bin_cache = channel_cache.at(i);
484 if( !bin_cache.hasStatUncert ) continue;
485 RooRealVar* gamma = bin_cache.gamma;
486 gamma->setVal(1.0);
487 }
488 std::vector< double > nu_b_stat_vec( channel_cache.size() );
489 for( unsigned int i = 0; i < channel_cache.size(); ++i ) {
490 BarlowCache& bin_cache = channel_cache.at(i);
491 if( !bin_cache.hasStatUncert ) continue;
492
493 RooAbsPdf* sum_pdf = (RooAbsPdf*) bin_cache.sumPdf;
494 RooArgSet* obsSet = bin_cache.observables;
495 double binVolume = bin_cache.binVolume;
496
497 bin_cache.SetBinCenter();
498 double nu_b_stat = sum_pdf->getVal(*obsSet)*sum_pdf->expectedEvents(*obsSet)*binVolume - nu_b_vec.at(i);
499 nu_b_stat_vec.at(i) = nu_b_stat;
500 }
501 //time_after_setVal=clock();
502
503 // Done with the first loops.
504 // Now evaluating the function
505
506 //clock_t time_before_eval, time_after_eval;
507
508 // Loop over the bins in the cache
509 //time_before_eval=clock();
510 for( unsigned int i = 0; i < channel_cache.size(); ++i ) {
511
512 BarlowCache& bin_cache = channel_cache.at(i);
513
514 if( !bin_cache.hasStatUncert ) {
515 //std::cout << "Bin: " << i << " of " << channel_cache.size()
516 // << " in channel: " << channel_name
517 // << " doesn't have stat uncertainties" << std::endl;
518 continue;
519 }
520
521 // Set the observable to the bin center
522 bin_cache.SetBinCenter();
523
524 // Get the cached objects
525 RooRealVar* gamma = bin_cache.gamma;
526 RooRealVar* tau = bin_cache.tau;
527 RooAbsReal* pois_mean = bin_cache.nom_pois_mean;
528 //RooAbsPdf* sum_pdf = (RooAbsPdf*) bin_cache.sumPdf;
529 //RooArgSet* obsSet = bin_cache.observables;
530 //double binVolume = bin_cache.binVolume;
531
532 // Get the values necessary for
533 // the analytic minimization
534 double nu_b = nu_b_vec.at(i);
535 double nu_b_stat = nu_b_stat_vec.at(i);
536
537 double tau_val = tau->getVal();
538 double nData = bin_cache.nData;
539 double m_val = pois_mean->getVal();
540
541 // Initialize the minimized value of gamma
542 double gamma_hat_hat = 1.0;
543
544 // Check that the quadratic term is > 0
545 if(nu_b_stat > 0.00000001) {
546
547 double A = nu_b_stat*nu_b_stat + tau_val*nu_b_stat;
548 double B = nu_b*tau_val + nu_b*nu_b_stat - nData*nu_b_stat - m_val*nu_b_stat;
549 double C = -1*m_val*nu_b;
550
551 double discrim = B*B-4*A*C;
552
553 if( discrim < 0 ) {
554 std::cout << "Warning: Discriminant (B*B - 4AC) < 0" << std::endl;
555 std::cout << "Warning: Taking B*B - 4*A*C == 0" << std::endl;
556 discrim=0;
557 //throw runtime_error("BarlowBeestonLL::evaluate() : B*B - 4AC < 0");
558 }
559 if( A <= 0 ) {
560 std::cout << "Warning: A <= 0" << std::endl;
561 throw runtime_error("BarlowBeestonLL::evaluate() : A < 0");
562 }
563
564 gamma_hat_hat = ( -1*B + TMath::Sqrt(discrim) ) / (2*A);
565 }
566
567 // If the quadratic term is 0, we simply
568 // use a linear equation
569 else {
570 gamma_hat_hat = m_val/tau_val;
571 }
572
573 // Check for NAN
574 if( TMath::IsNaN(gamma_hat_hat) ) {
575 std::cout << "ERROR: gamma hat hat is NAN" << std::endl;
576 throw runtime_error("BarlowBeestonLL::evaluate() : gamma hat hat is NAN");
577 }
578
579 if( gamma_hat_hat <= 0 ) {
580 std::cout << "WARNING: gamma hat hat <= 0. Setting to 0" << std::endl;
581 gamma_hat_hat = 0;
582 }
583
584 /*
585 std::cout << "n: " << bin_cache.nData << " "
586 << "nu_stat: " << nu_b_stat << " "
587 << "nu: " << nu_b << " "
588 << "tau: " << tau->getVal() << " "
589 << "m: " << pois_mean->getVal() << " "
590 << "A: " << A << " "
591 << "B: " << B << " "
592 << "C: " << C << " "
593 << "gamma hat hat: " << gamma_hat_hat
594 << std::endl;
595 */
596
597 gamma->setVal( gamma_hat_hat );
598
599 }
600
601 //time_after_eval=clock();
602
603 //float time_setVal = ((float) time_after_setVal - (float) time_before_setVal) / ((float) CLOCKS_PER_SEC);
604 //float time_eval = ((float) time_after_eval - (float) time_before_eval) / ((float) CLOCKS_PER_SEC);
605
606 /*
607 std::cout << "Barlow timing for channel: " << channel_name
608 << " SetVal: " << time_setVal
609 << " Eval: " << time_eval
610 << std::endl;
611 */
612 }
613
614
615 return _nll;
616
617}
618
619
620
621/*
622////////////////////////////////////////////////////////////////////////////////
623/// Check that parameters and likelihood value for 'best fit' are still valid. If not,
624/// because the best fit has never been calculated, or because constant parameters have
625/// changed value or parameters have changed const/float status, the minimum is recalculated
626
627void RooStats::HistFactory::RooBarlowBeestonLL::validateAbsMin() const
628{
629 // Check if constant status of any of the parameters have changed
630 if (_absMinValid) {
631 _piter->Reset() ;
632 RooAbsArg* par ;
633 while((par=(RooAbsArg*)_piter->Next())) {
634 if (_paramFixed[par->GetName()] != par->isConstant()) {
635 cxcoutI(Minimization) << "RooStats::HistFactory::RooBarlowBeestonLL::evaluate(" << GetName() << ") constant status of parameter " << par->GetName() << " has changed from "
636 << (_paramFixed[par->GetName()]?"fixed":"floating") << " to " << (par->isConstant()?"fixed":"floating")
637 << ", recalculating absolute minimum" << endl ;
638 _absMinValid = false ;
639 break ;
640 }
641 }
642 }
643
644
645 // If we don't have the absolute minimum w.r.t all observables, calculate that first
646 if (!_absMinValid) {
647
648 cxcoutI(Minimization) << "RooStats::HistFactory::RooBarlowBeestonLL::evaluate(" << GetName() << ") determining minimum likelihood for current configurations w.r.t all observable" << endl ;
649
650
651 // Save current values of non-marginalized parameters
652 RooArgSet* obsStart = (RooArgSet*) _obs.snapshot(false) ;
653
654 // Start from previous global minimum
655 if (_paramAbsMin.getSize()>0) {
656 const_cast<RooSetProxy&>(_par).assignValueOnly(_paramAbsMin) ;
657 }
658 if (_obsAbsMin.getSize()>0) {
659 const_cast<RooSetProxy&>(_obs).assignValueOnly(_obsAbsMin) ;
660 }
661
662 // Find minimum with all observables floating
663 const_cast<RooSetProxy&>(_obs).setAttribAll("Constant",false) ;
664 _minuit->migrad() ;
665
666 // Save value and remember
667 _absMin = _nll ;
668 _absMinValid = true ;
669
670 // Save parameter values at abs minimum as well
671 _paramAbsMin.removeAll() ;
672
673 // Only store non-constant parameters here!
674 RooArgSet* tmp = (RooArgSet*) _par.selectByAttrib("Constant",false) ;
675 _paramAbsMin.addClone(*tmp) ;
676 delete tmp ;
677
678 _obsAbsMin.addClone(_obs) ;
679
680 // Save constant status of all parameters
681 _piter->Reset() ;
682 RooAbsArg* par ;
683 while((par=(RooAbsArg*)_piter->Next())) {
684 _paramFixed[par->GetName()] = par->isConstant() ;
685 }
686
687 if (dologI(Minimization)) {
688 cxcoutI(Minimization) << "RooStats::HistFactory::RooBarlowBeestonLL::evaluate(" << GetName() << ") minimum found at (" ;
689
690 RooAbsReal* arg ;
691 bool first=true ;
692 _oiter->Reset() ;
693 while ((arg=(RooAbsReal*)_oiter->Next())) {
694 ccxcoutI(Minimization) << (first?"":", ") << arg->GetName() << "=" << arg->getVal() ;
695 first=false ;
696 }
697 ccxcoutI(Minimization) << ")" << endl ;
698 }
699
700 // Restore original parameter values
701 const_cast<RooSetProxy&>(_obs) = *obsStart ;
702 delete obsStart ;
703
704 }
705}
706*/
#define ClassImp(name)
Definition: Rtypes.h:375
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t target
char name[80]
Definition: TGX11.cxx:110
A class which maps the current values of a RooRealVar (or a set of RooRealVars) to one of a number of...
Definition: ParamHistFunc.h:24
const RooArgSet * get(Int_t masterIdx) const
Definition: ParamHistFunc.h:46
Int_t numBins() const
Definition: ParamHistFunc.h:36
RooAbsReal & getParameter() const
double binVolume() const
Definition: ParamHistFunc.h:49
bool isConstant() const
Check if the "Constant" attribute is set.
Definition: RooAbsArg.h:362
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...
Definition: RooAbsArg.cxx:541
virtual bool remove(const RooAbsArg &var, bool silent=false, bool matchByNameOnly=false)
Remove the specified argument from our list.
bool empty() const
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
void reserve(Storage_t::size_type count)
virtual double expectedEvents(const RooArgSet *nset) const
Return expected number of events to be used in calculation of extended likelihood.
Definition: RooAbsPdf.cxx:3195
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:62
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:91
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:56
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition: RooArgSet.h:179
RooCategory is an object to represent discrete states.
Definition: RooCategory.h:28
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:40
RooSimultaneous facilitates simultaneous fitting of multiple PDFs to subsets of a given dataset.
const RooAbsCategoryLValue & indexCat() const
RooAbsPdf * getPdf(const char *catName) const
Return the p.d.f associated with the given index category name.
Class RooBarlowBeestonLL implements the profile likelihood estimator for a given likelihood and set o...
double evaluate() const override
Optimized implementation of createProfile for profile likelihoods.
bool getParameters(const RooArgSet *depList, RooArgSet &outputSet, bool stripDisconnected=true) const override
Fills a list with leaf nodes in the arg tree starting with ourself as top node that don't match any o...
const char * GetName() const override
Returns name of object.
Definition: TNamed.h:47
bool getStatUncertaintyFromChannel(RooAbsPdf *channel, ParamHistFunc *&paramfunc, RooArgList *gammaList)
void FactorizeHistFactoryPdf(const RooArgSet &, RooAbsPdf &, RooArgList &, RooArgList &)
void getDataValuesForObservables(std::map< std::string, std::vector< double > > &ChannelBinDataMap, RooAbsData *data, RooAbsPdf *simPdf)
RooAbsPdf * getSumPdfFromChannel(RooAbsPdf *channel)
int getStatUncertaintyConstraintTerm(RooArgList *constraints, RooRealVar *gamma_stat, RooAbsReal *&pois_mean, RooRealVar *&tau)
static double B[]
static double A[]
static double C[]
double gamma(double x)
Bool_t IsNaN(Double_t x)
Definition: TMath.h:890
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition: TMath.h:660