#include <SimplerLikelihoodRatioTestStatExt.h>

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)
std::auto_ptr< RooArgSet >	paramsAlt_
std::auto_ptr< RooArgSet >	paramsNull_
	parameter sets to apply snapshots to
RooAbsPdf *	pdfAlt_
std::auto_ptr< RooAbsPdf >	pdfAltOwned_
RooArgSet	pdfCompAlt_
RooArgSet	pdfCompNull_
	snapshot with all branch nodes of the pdfs before factorization
RooArgList	pdfDepObs_
	nodes which depend directly on observables, on which one has to do redirectServers
RooAbsPdf *	pdfNull_
	pdfs (cloned, and with constraints factorized away)
std::auto_ptr< RooAbsPdf >	pdfNullOwned_
	owned copy of the pdfs after factorizing
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
RooSimultaneous *	simPdfNull_
	pdfNull, pdfAlt cast to sympdf (may be null) after factorization
RooArgSet	snapAlt_
RooArgSet	snapNull_
	snapshots of parameters for the two pdfs

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.

{
}

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, cmsDownloadME::cat, 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 [inline, virtual]

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;
        }
    }
}