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