ROOT  6.07/01
Reference Guide
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
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 //
13 // BEGIN_HTML
14 // Class RooBarlowBeestonLL implements the profile likelihood estimator for
15 // a given likelihood and set of parameters of interest. The value return by
16 // RooBarlowBeestonLL is the input likelihood nll minimized w.r.t all nuisance parameters
17 // (which are all parameters except for those listed in the constructor) minus
18 // the -log(L) of the best fit. Note that this function is slow to evaluate
19 // as a MIGRAD minimization step is executed for each function evaluation
20 // END_HTML
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 
38 #include "RooStats/RooStatsUtils.h"
39 #include "RooProdPdf.h"
40 #include "RooCategory.h"
41 #include "RooSimultaneous.h"
42 #include "RooArgList.h"
43 #include "RooAbsCategoryLValue.h"
44 
47 
48 using namespace std ;
49 
51 
52 
53 ////////////////////////////////////////////////////////////////////////////////
54 
55  RooStats::HistFactory::RooBarlowBeestonLL::RooBarlowBeestonLL() :
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 
399 const RooArgSet& RooStats::HistFactory::RooBarlowBeestonLL::bestFitParams() const
400 {
401  validateAbsMin() ;
402  return _paramAbsMin ;
403 }
404 
405 
406 ////////////////////////////////////////////////////////////////////////////////
407 
408 const 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 /*
423 RooAbsReal* RooStats::HistFactory::RooBarlowBeestonLL::createProfile(const RooArgSet& paramsOfInterest)
424 {
425  return nll().createProfile(paramsOfInterest) ;
426 }
427 */
428 
429 
430 /*
431 void 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 
666 void 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 
static double B[]
RooAbsCollection * snapshot(Bool_t deepCopy=kTRUE) const
Take a snap shot of current collection contents: An owning collection is returned containing clones o...
int getStatUncertaintyConstraintTerm(RooArgList *constraints, RooRealVar *gamma_stat, RooAbsReal *&pois_mean, RooRealVar *&tau)
bool getStatUncertaintyFromChannel(RooAbsPdf *channel, ParamHistFunc *&paramfunc, RooArgList *gammaList)
const RooAbsCategoryLValue & indexCat() const
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
const Bool_t kFALSE
Definition: Rtypes.h:92
RooAbsPdf * getPdf(const char *catName) const
Return the p.d.f associated with the given index category name.
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:2930
static double A[]
Iterator abstract base class.
Definition: TIterator.h:32
TText * tt
Definition: textangle.C:16
std::map< std::string, std::string >::const_iterator iter
Definition: TAlienJob.cxx:54
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...
TIterator * createIterator(Bool_t dir=kIterForward) const
RooCatType is an auxilary class for RooAbsCategory and defines a a single category state...
Definition: RooCatType.h:23
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
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:560
double binVolume() const
Definition: ParamHistFunc.h:54
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:37
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:203
double gamma(double x)
static double C[]
TPaveLabel title(3, 27.1, 15, 28.7,"ROOT Environment and Tools")
ClassImp(RooStats::HistFactory::RooBarlowBeestonLL) RooStats
bool verbose
void getDataValuesForObservables(std::map< std::string, std::vector< double > > &ChannelBinDataMap, RooAbsData *data, RooAbsPdf *simPdf)
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:51
RooCategory represents a fundamental (non-derived) discrete value object.
Definition: RooCategory.h:25
Bool_t isConstant() const
Definition: RooAbsArg.h:266
double Double_t
Definition: RtypesCore.h:55
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:53
#define name(a, b)
Definition: linkTestLib0.cpp:5
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.
const RooArgSet * get(Int_t masterIdx) const
Definition: ParamHistFunc.h:51
RooAbsPdf is the abstract interface for all probability density functions The class provides hybrid a...
Definition: RooAbsPdf.h:41
virtual TObject * Next()=0
Int_t IsNaN(Double_t x)
Definition: TMath.h:617
#define NULL
Definition: Rtypes.h:82
Int_t getSize() const
Double_t Sqrt(Double_t x)
Definition: TMath.h:464
RooAbsPdf * getSumPdfFromChannel(RooAbsPdf *channel)
const Bool_t kTRUE
Definition: Rtypes.h:91
RooSimultaneous facilitates simultaneous fitting of multiple PDFs to subsets of a given dataset...
RooRealVar & getParameter() const
Int_t numBins() const
Definition: ParamHistFunc.h:41
virtual const Text_t * GetName() const
Returns name of object.
Definition: RooCatType.h:45
void FactorizeHistFactoryPdf(const RooArgSet &, RooAbsPdf &, RooArgList &, RooArgList &)
TIterator * typeIterator() const
Return iterator over all defined states.
virtual Bool_t redirectServersHook(const RooAbsCollection &, Bool_t, Bool_t, Bool_t)
Double_t evaluate() const
Optimized implementation of createProfile for profile likelihoods.