54 RooAbsReal(
"RooBarlowBeestonLL",
"RooBarlowBeestonLL"),
58 _pdf(NULL), _data(NULL)
72 _nll(
"input",
"-log(L) function",this,nllIn),
75 _pdf(NULL), _data(NULL)
104 _nll(
"nll",this,other._nll),
107 _pdf(NULL), _data(NULL),
108 _paramFixed(other._paramFixed)
141 TIterator* iter = bin_center->createIterator() ;
157 std::cout <<
"Error: Must initialize data before initializing cache" << std::endl;
158 throw std::runtime_error(
"Uninitialized Data");
161 std::cout <<
"Error: Must initialize model pdf before initializing cache" << std::endl;
162 throw std::runtime_error(
"Uninitialized model pdf");
166 std::map< std::string, std::vector<double> > ChannelBinDataMap;
172 RooArgSet* obsSet = _pdf->getObservables(*_data);
175 if( obsTerms.
getSize() == 0 ) {
176 std::cout <<
"Error: Found no observable terms with pdf: " << _pdf->GetName()
177 <<
" using dataset: " << _data->GetName() << std::endl;
180 if( constraints.
getSize() == 0 ) {
181 std::cout <<
"Error: Found no constraint terms with pdf: " << _pdf->GetName()
182 <<
" using dataset: " << _data->GetName() << std::endl;
196 for (
const auto& nameIdx : *channelCat) {
201 std::string channel_name = channelPdf->
GetName();
207 if( ! hasStatUncert ) {
209 std::cout <<
"Channel: " << channel_name
210 <<
" doesn't have statistical uncertainties"
216 if(
verbose) std::cout <<
"Found ParamHistFunc: " << param_func->
GetName() << std::endl;
223 int num_bins = param_func->
numBins();
228 std::vector<BarlowCache> temp_cache( num_bins );
229 bool channel_has_stat_uncertainty=
false;
231 for(
Int_t bin_index = 0; bin_index < num_bins; ++bin_index ) {
239 if(
verbose) std::cout <<
"Ignoring constant gamma: " << gamma_stat->
GetName() << std::endl;
244 channel_has_stat_uncertainty=
true;
245 cache.
gamma = gamma_stat;
246 _statUncertParams.insert( gamma_stat->
GetName() );
263 if( !tau || !pois_mean ) {
264 std::cout <<
"Failed to find pois mean or tau parameter for " << gamma_stat->
GetName() << std::endl;
267 if(
verbose) std::cout <<
"Found pois mean and tau for parameter: " << gamma_stat->
GetName()
269 <<
" pois_mean: " << pois_mean->
GetName() <<
" " << pois_mean->
getVal()
278 if( sum_pdf == NULL ) {
279 std::cout <<
"Failed to find RooRealSumPdf in channel " << channel_name
280 <<
", therefor skipping this channel for analytic uncertainty minimization"
282 channel_has_stat_uncertainty=
false;
288 if( ChannelBinDataMap.find(channel_name) == ChannelBinDataMap.end() ) {
289 std::cout <<
"Error: channel with name: " << channel_name
290 <<
" not found in BinDataMap" << std::endl;
291 throw runtime_error(
"BinDataMap");
293 double nData = ChannelBinDataMap[channel_name].at(bin_index);
296 temp_cache.at(bin_index) = cache;
301 if( channel_has_stat_uncertainty ) {
302 std::cout <<
"Adding channel: " << channel_name
303 <<
" to the barlow cache" << std::endl;
304 _barlowCache[channel_name] = temp_cache;
365 std::string arg_name = arg->
GetName();
370 if( _statUncertParams.find(arg_name.c_str()) != _statUncertParams.end() ) {
384const RooArgSet& RooStats::HistFactory::RooBarlowBeestonLL::bestFitParams() const
387 return _paramAbsMin ;
393const RooArgSet& RooStats::HistFactory::RooBarlowBeestonLL::bestFitObs() const
465 std::map< std::string, std::vector< BarlowCache > >::iterator iter_cache;
466 for( iter_cache = _barlowCache.begin(); iter_cache != _barlowCache.end(); ++iter_cache ) {
468 std::string channel_name = (*iter_cache).first;
469 std::vector< BarlowCache >& channel_cache = (*iter_cache).second;
484 for(
unsigned int i = 0; i < channel_cache.size(); ++i ) {
490 std::vector< double > nu_b_vec( channel_cache.size() );
491 for(
unsigned int i = 0; i < channel_cache.size(); ++i ) {
501 nu_b_vec.at(i) = nu_b;
506 for(
unsigned int i = 0; i < channel_cache.size(); ++i ) {
512 std::vector< double > nu_b_stat_vec( channel_cache.size() );
513 for(
unsigned int i = 0; i < channel_cache.size(); ++i ) {
522 double nu_b_stat = sum_pdf->
getVal(*obsSet)*sum_pdf->
expectedEvents(*obsSet)*binVolume - nu_b_vec.at(i);
523 nu_b_stat_vec.at(i) = nu_b_stat;
534 for(
unsigned int i = 0; i < channel_cache.size(); ++i ) {
558 double nu_b = nu_b_vec.at(i);
559 double nu_b_stat = nu_b_stat_vec.at(i);
561 double tau_val = tau->
getVal();
562 double nData = bin_cache.
nData;
563 double m_val = pois_mean->
getVal();
566 double gamma_hat_hat = 1.0;
569 if(nu_b_stat > 0.00000001) {
571 double A = nu_b_stat*nu_b_stat + tau_val*nu_b_stat;
572 double B = nu_b*tau_val + nu_b*nu_b_stat - nData*nu_b_stat - m_val*nu_b_stat;
573 double C = -1*m_val*nu_b;
575 double discrim =
B*
B-4*
A*
C;
578 std::cout <<
"Warning: Discriminant (B*B - 4AC) < 0" << std::endl;
579 std::cout <<
"Warning: Taking B*B - 4*A*C == 0" << std::endl;
584 std::cout <<
"Warning: A <= 0" << std::endl;
585 throw runtime_error(
"BarlowBeestonLL::evaluate() : A < 0");
594 gamma_hat_hat = m_val/tau_val;
599 std::cout <<
"ERROR: gamma hat hat is NAN" << std::endl;
600 throw runtime_error(
"BarlowBeestonLL::evaluate() : gamma hat hat is NAN");
603 if( gamma_hat_hat <= 0 ) {
604 std::cout <<
"WARNING: gamma hat hat <= 0. Setting to 0" << std::endl;
621 gamma->setVal( gamma_hat_hat );
651void RooStats::HistFactory::RooBarlowBeestonLL::validateAbsMin() const
653 // Check if constant status of any of the parameters have changed
657 while((par=(RooAbsArg*)_piter->Next())) {
658 if (_paramFixed[par->GetName()] != par->isConstant()) {
659 cxcoutI(Minimization) << "RooStats::HistFactory::RooBarlowBeestonLL::evaluate(" << GetName() << ") constant status of parameter " << par->GetName() << " has changed from "
660 << (_paramFixed[par->GetName()]?"fixed":"floating") << " to " << (par->isConstant()?"fixed":"floating")
661 << ", recalculating absolute minimum" << endl ;
662 _absMinValid = kFALSE ;
669 // If we don't have the absolute minimum w.r.t all observables, calculate that first
672 cxcoutI(Minimization) << "RooStats::HistFactory::RooBarlowBeestonLL::evaluate(" << GetName() << ") determining minimum likelihood for current configurations w.r.t all observable" << endl ;
675 // Save current values of non-marginalized parameters
676 RooArgSet* obsStart = (RooArgSet*) _obs.snapshot(kFALSE) ;
678 // Start from previous global minimum
679 if (_paramAbsMin.getSize()>0) {
680 const_cast<RooSetProxy&>(_par).assignValueOnly(_paramAbsMin) ;
682 if (_obsAbsMin.getSize()>0) {
683 const_cast<RooSetProxy&>(_obs).assignValueOnly(_obsAbsMin) ;
686 // Find minimum with all observables floating
687 const_cast<RooSetProxy&>(_obs).setAttribAll("Constant",kFALSE) ;
690 // Save value and remember
692 _absMinValid = kTRUE ;
694 // Save parameter values at abs minimum as well
695 _paramAbsMin.removeAll() ;
697 // Only store non-constant parameters here!
698 RooArgSet* tmp = (RooArgSet*) _par.selectByAttrib("Constant",kFALSE) ;
699 _paramAbsMin.addClone(*tmp) ;
702 _obsAbsMin.addClone(_obs) ;
704 // Save constant status of all parameters
707 while((par=(RooAbsArg*)_piter->Next())) {
708 _paramFixed[par->GetName()] = par->isConstant() ;
711 if (dologI(Minimization)) {
712 cxcoutI(Minimization) << "RooStats::HistFactory::RooBarlowBeestonLL::evaluate(" << GetName() << ") minimum found at (" ;
717 while ((arg=(RooAbsReal*)_oiter->Next())) {
718 ccxcoutI(Minimization) << (first?"":", ") << arg->GetName() << "=" << arg->getVal() ;
721 ccxcoutI(Minimization) << ")" << endl ;
724 // Restore original parameter values
725 const_cast<RooSetProxy&>(_obs) = *obsStart ;
A class which maps the current values of a RooRealVar (or a set of RooRealVars) to one of a number of...
const RooArgSet * get(Int_t masterIdx) const
RooRealVar & getParameter() const
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...
Bool_t isConstant() const
Check if the "Constant" attribute is set.
RooAbsCollection is an abstract container object that can hold multiple RooAbsArg objects.
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
TIterator-style iteration over contained elements.
virtual Double_t expectedEvents(const RooArgSet *nset) const
Return expected number of events from this p.d.f for use in extended likelihood calculations.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
RooArgList is a container object that can hold multiple RooAbsArg objects.
RooArgSet is a container object that can hold multiple RooAbsArg objects.
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
RooCategory is an object to represent discrete states.
RooRealVar represents a variable that can be changed from the outside.
virtual void setVal(Double_t value)
Set value of variable to 'value'.
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.
RooAbsReal * nom_pois_mean
void SetBinCenter() const
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.
void initializeBarlowCache()
virtual ~RooBarlowBeestonLL()
Destructor.
Iterator abstract base class.
virtual TObject * Next()=0
virtual const char * GetName() const
Returns name of object.
bool getStatUncertaintyFromChannel(RooAbsPdf *channel, ParamHistFunc *¶mfunc, 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)
Double_t Sqrt(Double_t x)