#include <SimplerLikelihoodRatioTestStatExt.h>

Inheritance diagram for SimplerLikelihoodRatioTestStatOpt:

Public Member Functions
virtual Double_t	Evaluate (RooAbsData &data, RooArgSet &nullPOI)

virtual const TString	GetVarName () const

	SimplerLikelihoodRatioTestStatOpt (const RooArgSet &obs, RooAbsPdf &pdfNull, RooAbsPdf &pdfAlt, const RooArgSet &paramsNull=RooArgSet(), const RooArgSet &paramsAlt=RooArgSet(), bool factorize=true)

virtual	~SimplerLikelihoodRatioTestStatOpt ()

Private Member Functions
double	evalSimNLL (RooAbsData &data, RooSimultaneous pdf, std::vector< RooAbsPdf > &components)

double	evalSimpleNLL (RooAbsData &data, RooAbsPdf *pdf)

void	unrollSimPdf (RooSimultaneous pdf, std::vector< RooAbsPdf > &out)

Private Attributes
const RooArgSet *	obs_
	observables (global argset) More...

std::auto_ptr< RooArgSet >	paramsAlt_

std::auto_ptr< RooArgSet >	paramsNull_
	parameter sets to apply snapshots to More...

RooAbsPdf *	pdfAlt_

std::auto_ptr< RooAbsPdf >	pdfAltOwned_

RooArgSet	pdfCompAlt_

RooArgSet	pdfCompNull_
	snapshot with all branch nodes of the pdfs before factorization More...

RooArgList	pdfDepObs_
	nodes which depend directly on observables, on which one has to do redirectServers More...

RooAbsPdf *	pdfNull_
	pdfs (cloned, and with constraints factorized away) More...

std::auto_ptr< RooAbsPdf >	pdfNullOwned_
	owned copy of the pdfs after factorizing More...

RooSimultaneous *	simPdfAlt_

std::vector< RooAbsPdf * >	simPdfComponentsAlt_

std::vector< RooAbsPdf * >	simPdfComponentsNull_
	components of the sim pdfs after factorization, for each bin in sim. category. can contain nulls More...

RooSimultaneous *	simPdfNull_
	pdfNull, pdfAlt cast to sympdf (may be null) after factorization More...

RooArgSet	snapAlt_

RooArgSet	snapNull_
	snapshots of parameters for the two pdfs More...

Detailed Description

BEGIN_HTML

SimplerLikelihoodRatioTestStatOpt: an optimized implementaion of the simple likelihood ratio test statistics Q = - ln( L(data|null) / L(data|alt) )

Optimization is obtained through the following methods:

Constraints not depending on observables are factorized away at construction time.
The nodes that directly depend on the observables are identified at construction time, and so only those nodes have their servers redirected when changing dataset.
Simultaneous pdfs are split at construction time; when evaluating the likelihood, no lookups of components by name are performed.
The amount of objects created and destroyed for each evaluation is reduced to a minimum.
Only one loop is performed on the dataset (since only one is necessary), and the data is not copied.

Author: Giovanni Petrucciani (UCSD/CMS/CERN), May 2011

END_HTML

Definition at line 34 of file SimplerLikelihoodRatioTestStatExt.h.

Constructor & Destructor Documentation

SimplerLikelihoodRatioTestStatOpt::SimplerLikelihoodRatioTestStatOpt	(	const RooArgSet &	obs,
		RooAbsPdf &	pdfNull,
		RooAbsPdf &	pdfAlt,
		const RooArgSet &	paramsNull = `RooArgSet()`,
		const RooArgSet &	paramsAlt = `RooArgSet()`,
		bool	factorize = `true`
	)

Create a SimplerLikelihoodRatioTestStatOpt. obs = the observables on which the two pdfs depend pdfNull, pdfAlt = pdfs of the two models paramsNull, paramsAlt = values of the parameters that should be set before evaluating the pdfs factorize: if set to true, the constraint terms not depending on the observables will be removed

Definition at line 4 of file SimplerLikelihoodRatioTestStatExt.cc.

References createBeamHaloJobs::constraints, utils::factorizePdf(), utils::fullClonePdf(), utils::getClients(), obs_, pdfAlt_, pdfAltOwned_, pdfCompAlt_, pdfCompNull_, pdfDepObs_, pdfNull_, pdfNullOwned_, simPdfAlt_, simPdfComponentsAlt_, simPdfComponentsNull_, simPdfNull_, snapAlt_, snapNull_, and unrollSimPdf().

                         :
     obs_(&obs)
 {
     // take a snapshot of the pdf so we can modify it at will
     RooAbsPdf *cloneNull = utils::fullClonePdf(&pdfNull, pdfCompNull_);
     RooAbsPdf *cloneAlt  = utils::fullClonePdf(&pdfAlt, pdfCompAlt_);
 
     // factorize away constraints
     RooArgList constraints;
     pdfNull_ = factorize ? utils::factorizePdf(obs, *cloneNull, constraints) : cloneNull;
     if (pdfNull_ == 0) throw std::invalid_argument("SimplerLikelihoodRatioTestStatOpt:: pdf does not depend on observables");
     if (pdfNull_ != cloneNull) pdfNullOwned_.reset(pdfNull_); // if new, then take ownership of it
     if (cloneNull == cloneAlt) {
         pdfAlt_ = pdfNull_;
     } else {
         pdfAlt_  = factorize ? utils::factorizePdf(obs, *cloneAlt,  constraints) : cloneAlt;
         if (pdfAlt_ == 0) throw std::invalid_argument("SimplerLikelihoodRatioTestStatOpt:: pdf does not depend on observables");
         if (pdfAlt_ != cloneAlt)  pdfAltOwned_.reset(pdfAlt_); // if new, then take ownership of it
     }
 
     // take snapshot of parameters
     snapNull_.addClone(paramsNull);
     snapAlt_.addClone(paramsAlt);
 
     // if the pdf is a RooSimultaneous, unroll it
     simPdfNull_ =  dynamic_cast<RooSimultaneous *>(pdfNull_);
     simPdfAlt_ =  dynamic_cast<RooSimultaneous *>(pdfAlt_);
     if (simPdfNull_) unrollSimPdf(simPdfNull_, simPdfComponentsNull_);
     if (simPdfAlt_) unrollSimPdf(simPdfAlt_, simPdfComponentsAlt_);
 
     // optimize caching
     pdfNull_->optimizeCacheMode(*obs_);
     pdfAlt_->optimizeCacheMode(*obs_);
 
     // find nodes that directly depend on observables
     utils::getClients(*obs_, pdfCompNull_, pdfDepObs_);
     utils::getClients(*obs_, pdfCompAlt_,  pdfDepObs_);
 }

SimplerLikelihoodRatioTestStatOpt::~SimplerLikelihoodRatioTestStatOpt ( )

virtual

Definition at line 47 of file SimplerLikelihoodRatioTestStatExt.cc.

48 {

49 }

Member Function Documentation

double SimplerLikelihoodRatioTestStatOpt::evalSimNLL	(	RooAbsData &	data,
		RooSimultaneous *	pdf,
		std::vector< RooAbsPdf * > &	components
	)

private

Definition at line 106 of file SimplerLikelihoodRatioTestStatExt.cc.

References newFWLiteAna::bin, i, n, obs_, and w().

Referenced by Evaluate().

                                                                                                                                 {
     data.setDirtyProp(false);
 
     double sum = 0.0;
     int i, n = data.numEntries(); 
 
     // must fetch the category, from the dataset first, if it's not empty
     RooAbsCategoryLValue *cat = 0;
     if (n) {
         const RooArgSet *entry = data.get(0);
         cat = dynamic_cast<RooAbsCategoryLValue *>(entry->find(pdf->indexCat().GetName()));
         assert(cat != 0 && "Didn't find category in dataset");
     }
 
     // now loop on the dataset, and dispatch the request to the appropriate pdf
     std::vector<double> sumw(components.size(), 0);
     for (i = 0; i < n; ++i) {
         data.get(i); 
         double w = data.weight(); if (w == 0) continue;
         int bin = cat->getBin();
         assert(bin < int(components.size()) && "Bin outside range");
         if (components[bin] == 0) continue;
         sum  += -w*components[bin]->getLogVal(obs_);
         sumw[bin] +=  w;
     }
 
     // then compute extended term
     for (i = 0, n = components.size(); i < n; ++i) {
         if (components[i]) sum += components[i]->extendedTerm(UInt_t(sumw[i]), obs_);
     }
     return sum;
 }

double SimplerLikelihoodRatioTestStatOpt::evalSimpleNLL	(	RooAbsData &	data,
		RooAbsPdf *	pdf
	)

private

Definition at line 139 of file SimplerLikelihoodRatioTestStatExt.cc.

References i, n, obs_, and w().

Referenced by Evaluate().

                                                                                          {
     data.setDirtyProp(false);
     double sum = 0.0, sumw = 0.0;
     int i, n = data.numEntries(); 
     for (i = 0; i < n; ++i) {
         data.get(i); 
         double w = data.weight(); if (w == 0) continue;
         sum  += -w*pdf->getLogVal(obs_);
         sumw +=  w;
     }
     sum += pdf->extendedTerm(UInt_t(sumw), obs_);
     return sum;
 }

Double_t SimplerLikelihoodRatioTestStatOpt::Evaluate	(	RooAbsData &	data,
		RooArgSet &	nullPOI
	)

virtual

Definition at line 52 of file SimplerLikelihoodRatioTestStatExt.cc.

References a, evalSimNLL(), evalSimpleNLL(), obs_, paramsAlt_, paramsNull_, pdfAlt_, pdfDepObs_, pdfNull_, simPdfAlt_, simPdfComponentsAlt_, simPdfComponentsNull_, simPdfNull_, snapAlt_, and snapNull_.

 {
     // get parameters, if not there already
     if (paramsNull_.get() == 0) paramsNull_.reset(pdfNull_->getParameters(data));
     if (paramsAlt_.get() == 0)  paramsAlt_.reset(pdfAlt_->getParameters(data));
 
     // if the dataset is not empty, redirect pdf nodes to the dataset
     std::auto_ptr<TIterator> iterDepObs(pdfDepObs_.createIterator());
     bool nonEmpty = data.numEntries() > 0;
     if (nonEmpty) {
         const RooArgSet *entry = data.get(0);
         for (RooAbsArg *a = (RooAbsArg *) iterDepObs->Next(); a != 0; a = (RooAbsArg *) iterDepObs->Next()) {
             a->redirectServers(*entry);    
         }
     }
 
     // evaluate null pdf
     *paramsNull_ = snapNull_;
     *paramsNull_ = nullPOI;
     double nullNLL = simPdfNull_ ? evalSimNLL(data, simPdfNull_, simPdfComponentsNull_) : evalSimpleNLL(data, pdfNull_);
 
     // evaluate alt pdf
     *paramsAlt_ = snapAlt_;
     double altNLL = simPdfAlt_ ? evalSimNLL(data, simPdfAlt_, simPdfComponentsAlt_) : evalSimpleNLL(data, pdfAlt_);
 
     // put back links in pdf nodes, otherwise if the dataset goes out of scope they have dangling pointers
     if (nonEmpty) {
         iterDepObs->Reset();
         for (RooAbsArg *a = (RooAbsArg *) iterDepObs->Next(); a != 0; a = (RooAbsArg *) iterDepObs->Next()) {
             a->redirectServers(*obs_);    
         }
     }
 
     return nullNLL-altNLL;
 }

virtual const TString SimplerLikelihoodRatioTestStatOpt::GetVarName ( ) const

inlinevirtual

Definition at line 47 of file SimplerLikelihoodRatioTestStatExt.h.

References pdfAlt_, and pdfNull_.

                                                  {
             return TString::Format("-log(%s/%s)", pdfNull_->GetName(), pdfAlt_->GetName()); 
         }

void SimplerLikelihoodRatioTestStatOpt::unrollSimPdf	(	RooSimultaneous *	pdf,
		std::vector< RooAbsPdf * > &	out
	)

private

Definition at line 88 of file SimplerLikelihoodRatioTestStatExt.cc.

References NULL, and obs_.

Referenced by SimplerLikelihoodRatioTestStatOpt().

                                                                                                          {
     // get a clone of the pdf category, so that I can use it to enumerate the pdf states
     std::auto_ptr<RooAbsCategoryLValue> catClone((RooAbsCategoryLValue*) simpdf->indexCat().Clone());
     out.resize(catClone->numBins(NULL), 0);
     //std::cout << "Pdf " << pdf->GetName() <<" is a SimPdf over category " << catClone->GetName() << ", with " << out.size() << " bins" << std::endl;
 
     // loop on the category state and fetch the pdfs
     for (int ib = 0, nb = out.size(); ib < nb; ++ib) {
         catClone->setBin(ib);
         RooAbsPdf *pdfi = simpdf->getPdf(catClone->getLabel());
         pdfi->optimizeCacheMode(*obs_);
         if (pdfi != 0) {
             //std::cout << "   bin " << ib << " (label " << catClone->getLabel() << ") has pdf " << pdfi->GetName() << " of type " << pdfi->ClassName() << std::endl;
             out[ib] = pdfi;
         }
     }
 }