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