#include <ProfiledLikelihoodRatioTestStatExt.h>

Public Types
enum	OneSidedness { twoSidedDef = 0, oneSidedDef = 1, signFlipDef = 2 }
Public Member Functions
virtual Double_t	Evaluate (RooAbsData &data, RooArgSet &nullPOI)
virtual std::vector< Double_t >	Evaluate (RooAbsData &data, RooArgSet &nullPOI, const std::vector< Double_t > &rVals)
virtual const TString	GetVarName () const
	ProfiledLikelihoodTestStatOpt (const RooArgSet &observables, RooAbsPdf &pdf, const RooArgSet *nuisances, const RooArgSet &params, const RooArgSet &poi, const RooArgList &gobsParams, const RooArgList &gobs, int verbosity=0, OneSidedness oneSided=oneSidedDef)
void	SetOneSided (OneSidedness oneSided)
void	setPrintLevel (Int_t level)
Private Member Functions
bool	createNLL (RooAbsPdf &pdf, RooAbsData &data)
double	minNLL (bool constrained, RooRealVar *r=0)
Private Attributes
RooArgList	gobs_
RooArgList	gobsParams_
std::auto_ptr< RooAbsReal >	nll_
RooArgSet	nuisances_
OneSidedness	oneSided_
std::auto_ptr< RooArgSet >	params_
RooAbsPdf *	pdf_
RooArgSet	poi_
RooArgSet	poiParams_
RooArgSet	snap_
Int_t	verbosity_

Detailed Description

Definition at line 49 of file ProfiledLikelihoodRatioTestStatExt.h.

Member Enumeration Documentation

enum ProfiledLikelihoodTestStatOpt::OneSidedness

Enumerator:

twoSidedDef
oneSidedDef
signFlipDef

Definition at line 51 of file ProfiledLikelihoodRatioTestStatExt.h.

{ twoSidedDef = 0, oneSidedDef = 1, signFlipDef = 2 };

Constructor & Destructor Documentation

ProfiledLikelihoodTestStatOpt::ProfiledLikelihoodTestStatOpt	(	const RooArgSet &	observables,
		RooAbsPdf &	pdf,
		const RooArgSet *	nuisances,
		const RooArgSet &	params,
		const RooArgSet &	poi,
		const RooArgList &	gobsParams,
		const RooArgList &	gobs,
		int	verbosity = `0`,
		OneSidedness	oneSided = `oneSidedDef`
	)

Definition at line 117 of file ProfiledLikelihoodRatioTestStatExt.cc.

References a, dtNoiseDBValidation_cfg::cerr, gather_cfg::cout, DBG, nuisances_, params_, pdf_, poi_, poiParams_, createTree::pp, and snap_.

:
    pdf_(&pdf),
    gobsParams_(gobsParams),
    gobs_(gobs),
    verbosity_(verbosity),
    oneSided_(oneSided)
{
    DBG(DBG_PLTestStat_main, (std::cout << "Created for " << pdf.GetName() << "." << std::endl))

    params.snapshot(snap_,false);
    ((RooRealVar*)snap_.find(params.first()->GetName()))->setConstant(false);
    if (nuisances) { nuisances_.add(*nuisances); snap_.addClone(*nuisances, /*silent=*/true); }
    params_.reset(pdf_->getParameters(observables));
    DBG(DBG_PLTestStat_ctor, (std::cout << "Observables: " << std::endl)) DBG(DBG_PLTestStat_ctor, (observables.Print("V")))
    DBG(DBG_PLTestStat_ctor, (std::cout << "All params: " << std::endl))  DBG(DBG_PLTestStat_ctor, (params_->Print("V")))
    DBG(DBG_PLTestStat_ctor, (std::cout << "Snapshot: " << std::endl))    DBG(DBG_PLTestStat_ctor, (snap_.Print("V")))
    DBG(DBG_PLTestStat_ctor, (std::cout << "POI: " << std::endl))         DBG(DBG_PLTestStat_ctor, (poi_.Print("V")))
    RooLinkedListIter it = poi.iterator();
    for (RooAbsArg *a = (RooAbsArg*) it.Next(); a != 0; a = (RooAbsArg*) it.Next()) {
        // search for this poi in the parameters and in the snapshot
        RooAbsArg *ps = snap_.find(a->GetName());   
        RooAbsArg *pp = params_->find(a->GetName());
        if (pp == 0) { std::cerr << "WARNING: NLL does not depend on POI " << a->GetName() << ", cannot profile" << std::endl; continue; }
        if (ps == 0) { std::cerr << "WARNING: no snapshot for POI " << a->GetName() << ", cannot profile"  << std::endl; continue; }
        poi_.add(*ps);
        poiParams_.add(*pp);
    }
}

Member Function Documentation

bool ProfiledLikelihoodTestStatOpt::createNLL	(	RooAbsPdf &	pdf,
		RooAbsData &	data
	)		`[private]`

Definition at line 362 of file ProfiledLikelihoodRatioTestStatExt.cc.

References nll_, and nuisances_.

Referenced by Evaluate().

{
    if (typeid(pdf) == typeid(RooSimultaneousOpt)) {
        if (nll_.get() == 0) nll_.reset(pdf.createNLL(data, RooFit::Constrain(nuisances_)));
        else ((cacheutils::CachingSimNLL&)(*nll_)).setData(data);
        return true;
    } else {
        nll_.reset(pdf.createNLL(data, RooFit::Constrain(nuisances_)));
        return false;
    }
}

std::vector< Double_t > ProfiledLikelihoodTestStatOpt::Evaluate	(	RooAbsData &	data,
		RooArgSet &	nullPOI,
		const std::vector< Double_t > &	rVals
	)		`[virtual]`

Definition at line 262 of file ProfiledLikelihoodRatioTestStatExt.cc.

References gather_cfg::cout, createNLL(), DBG, EPS, reco::get(), gobs_, gobsParams_, i, minNLL(), nll_, nuisances_, oneSided_, oneSidedDef, params_, pdf_, poi_, alignCSCRings::r, run_regression::ret, signFlipDef, and snap_.

{
    static bool do_debug = runtimedef::get("DEBUG_PLTSO");
    std::vector<Double_t> ret(rVals.size(), 0);

    DBG(DBG_PLTestStat_main, (std::cout << "Being evaluated on " << data.GetName() << std::endl))

    // Take snapshot of initial state, to restore it at the end 
    RooArgSet initialState; params_->snapshot(initialState);

    DBG(DBG_PLTestStat_pars, std::cout << "Being evaluated on " << data.GetName() << ": params before snapshot are " << std::endl)
    DBG(DBG_PLTestStat_pars, params_->Print("V"))
    // Initialize parameters
    *params_ = snap_;

    //set value of "globalConstrained" nuisances to the value of the corresponding global observable and set them constant
    for (int i=0; i<gobsParams_.getSize(); ++i) {
      RooRealVar *nuis = (RooRealVar*)gobsParams_.at(i);
      RooRealVar *gobs = (RooRealVar*)gobs_.at(i);
      nuis->setVal(gobs->getVal());
      nuis->setConstant();
    }
        
    DBG(DBG_PLTestStat_pars, std::cout << "Being evaluated on " << data.GetName() << ": params after snapshot are " << std::endl)
    DBG(DBG_PLTestStat_pars, params_->Print("V"))
    // Initialize signal strength
    RooRealVar *rIn = (RooRealVar *) poi_.first();
    RooRealVar *r   = (RooRealVar *) params_->find(rIn->GetName());
    bool canKeepNLL = createNLL(*pdf_, data);

    double initialR = rVals.back();

    // Perform unconstrained minimization (denominator)
    r->setMin(0); if (initialR == 0 || (oneSided_ != oneSidedDef)) r->removeMax(); else r->setMax(1.1*initialR); 
    r->setVal(initialR == 0 ? 0.5 : 0.5*initialR); //best guess
    r->setConstant(false);
    DBG(DBG_PLTestStat_pars, (std::cout << "r In: ")) DBG(DBG_PLTestStat_pars, (rIn->Print(""))) DBG(DBG_PLTestStat_pars, std::cout << std::endl)
    DBG(DBG_PLTestStat_pars, std::cout << "r before the fit: ") DBG(DBG_PLTestStat_pars, r->Print("")) DBG(DBG_PLTestStat_pars, std::cout << std::endl)

    if (do_debug) std::cout << "PERFORMING UNCONSTRAINED FIT " << r->GetName() << " [ " << r->getMin() << " - " << r->getVal() << " - " << r->getMax() << " ] "<< std::endl;
    double nullNLL = minNLL(/*constrained=*/false, r);
    double bestFitR = r->getVal();
    // Take snapshot of initial state, to restore it at the end 
    RooArgSet bestFitState; params_->snapshot(bestFitState);


    DBG(DBG_PLTestStat_pars, (std::cout << "r after the fit: ")) DBG(DBG_PLTestStat_pars, (r->Print(""))) DBG(DBG_PLTestStat_pars, std::cout << std::endl)
    DBG(DBG_PLTestStat_pars, std::cout << "Was evaluated on " << data.GetName() << ": params before snapshot are " << std::endl)
    DBG(DBG_PLTestStat_pars, params_->Print("V"))

    double EPS = 0.25*ROOT::Math::MinimizerOptions::DefaultTolerance();
    for (int iR = 0, nR = rVals.size(); iR < nR; ++iR) {
        if (fabs(ret[iR]) > 10*EPS) continue; // don't bother re-update points which were too far from zero anyway.
        *params_ = bestFitState;
        initialR = rVals[iR];
        // Prepare for constrained minimization (numerator)
        r->setVal(initialR); 
        r->setConstant(true);
        double thisNLL = nullNLL, sign = +1.0;
        if (initialR == 0 || oneSided_ != oneSidedDef || bestFitR < initialR) { 
            // must do constrained fit (if there's something to fit besides XS)
            //std::cout << "PERFORMING CONSTRAINED FIT " << r->GetName() << " == " << r->getVal() << std::endl;
            thisNLL = (nuisances_.getSize() > 0 ? minNLL(/*constrained=*/true, r) : nll_->getVal());
            if (thisNLL - nullNLL < 0 && thisNLL - nullNLL >= -EPS) {
                thisNLL = nullNLL;
            } else if (thisNLL - nullNLL < 0) {
                DBG(DBG_PLTestStat_main, (printf("  --> constrained fit is better... will repeat unconstrained fit\n")))
                r->setConstant(false);
                r->setVal(bestFitR);
                double oldNullNLL = nullNLL;
                nullNLL = minNLL(/*constrained=*/false, r);
                bestFitR = r->getVal();
                bestFitState.removeAll(); params_->snapshot(bestFitState);
                for (int iR2 = 0; iR2 < iR; ++iR2) {
                    ret[iR2] -= (nullNLL - oldNullNLL); // fixup already computed test statistics
                }
                iR = -1; continue; // restart over again, refitting those close to zero :-(
            }
            if (bestFitR > initialR && oneSided_ == signFlipDef) {
                DBG(DBG_PLTestStat_main, (printf("   fitted signal %7.4f > %7.4f, test statistics will be negative.\n", bestFitR, initialR)))
                sign = -1.0;
            }
        } else {
            DBG(DBG_PLTestStat_main, (printf("   signal fit %7.4f > %7.4f: don't need to compute numerator\n", bestFitR, initialR)))
        }

        ret[iR] = sign * (thisNLL-nullNLL);
        DBG(DBG_PLTestStat_main, (printf("\nNLLs for %7.4f:  num % 10.4f, den % 10.4f (signal %7.4f), test stat % 10.4f\n", initialR, thisNLL, nullNLL, bestFitR, ret[iR])))
        if (do_debug) { 
            printf("   Q(%.4f) = %.4f (best fit signal %.4f), from num %.4f, den %.4f\n", rVals[iR], ret[iR], bestFitR, thisNLL, nullNLL); fflush(stdout);
        }
    }
    //Restore initial state, to avoid issues with ToyMCSampler
    *params_ = initialState;

    if (!canKeepNLL) nll_.reset();

    //return std::min(thisNLL-nullNLL, 0.);
    return ret;
}

Double_t ProfiledLikelihoodTestStatOpt::Evaluate	(	RooAbsData &	data,
		RooArgSet &	nullPOI
	)		`[virtual]`

Definition at line 157 of file ProfiledLikelihoodRatioTestStatExt.cc.

References gather_cfg::cout, createNLL(), DBG, reco::get(), gobs_, gobsParams_, i, minNLL(), nll_, nuisances_, oneSided_, oneSidedDef, params_, pdf_, poi_, poiParams_, alignCSCRings::r, utils::setAllConstant(), signFlipDef, snap_, and swap().

Referenced by ProfileLikelihood::runSignificance().

{
    bool do_debug = runtimedef::get("DEBUG_PLTSO");
    //static bool do_rescan = runtimedef::get("PLTSO_FULL_RESCAN");
    DBG(DBG_PLTestStat_main, (std::cout << "Being evaluated on " << data.GetName() << std::endl))

    // Take snapshot of initial state, to restore it at the end 
    RooArgSet initialState; params_->snapshot(initialState);

    DBG(DBG_PLTestStat_pars, std::cout << "Being evaluated on " << data.GetName() << ": params before snapshot are " << std::endl)
    DBG(DBG_PLTestStat_pars, params_->Print("V"))
    // Initialize parameters
    *params_ = snap_;

    //set value of "globalConstrained" nuisances to the value of the corresponding global observable and set them constant
    for (int i=0; i<gobsParams_.getSize(); ++i) {
      RooRealVar *nuis = (RooRealVar*)gobsParams_.at(i);
      RooRealVar *gobs = (RooRealVar*)gobs_.at(i);
      nuis->setVal(gobs->getVal());
      nuis->setConstant();
    }
        
    DBG(DBG_PLTestStat_pars, std::cout << "Being evaluated on " << data.GetName() << ": params after snapshot are " << std::endl)
    DBG(DBG_PLTestStat_pars, params_->Print("V"))
    // Initialize signal strength (or the first parameter)
    RooRealVar *rIn = (RooRealVar *) poi_.first();
    RooRealVar *r   = (RooRealVar *) params_->find(rIn->GetName());
    bool canKeepNLL = createNLL(*pdf_, data);

    double initialR = rIn->getVal();

    // Perform unconstrained minimization (denominator)
    if (poi_.getSize() == 1) {
        r->setMin(0); if (initialR == 0 || (oneSided_ != oneSidedDef)) r->removeMax(); else r->setMax(1.1*initialR); 
        r->setVal(initialR == 0 ? 0.5 : 0.5*initialR); //best guess
        r->setConstant(false);
    } else {
        utils::setAllConstant(poiParams_,false);
    }
    DBG(DBG_PLTestStat_pars, (std::cout << "r In: ")) DBG(DBG_PLTestStat_pars, (rIn->Print(""))) DBG(DBG_PLTestStat_pars, std::cout << std::endl)
    DBG(DBG_PLTestStat_pars, std::cout << "r before the fit: ") DBG(DBG_PLTestStat_pars, r->Print("")) DBG(DBG_PLTestStat_pars, std::cout << std::endl)

    //std::cout << "PERFORMING UNCONSTRAINED FIT " << r->GetName() << " [ " << r->getMin() << " - " << r->getMax() << " ] "<< std::endl;
    double nullNLL = minNLL(/*constrained=*/false, r);
    double bestFitR = r->getVal();

    DBG(DBG_PLTestStat_pars, (std::cout << "r after the fit: ")) DBG(DBG_PLTestStat_pars, (r->Print(""))) DBG(DBG_PLTestStat_pars, std::cout << std::endl)
    DBG(DBG_PLTestStat_pars, std::cout << "Was evaluated on " << data.GetName() << ": params before snapshot are " << std::endl)
    DBG(DBG_PLTestStat_pars, params_->Print("V"))

    // Prepare for constrained minimization (numerator)
    if (poi_.getSize() == 1) {
        r->setVal(initialR); 
        r->setConstant(true);
    } else {
        poiParams_.assignValueOnly(poi_);
        utils::setAllConstant(poiParams_,true);
    }
    double thisNLL = nullNLL;
    if (initialR == 0 || oneSided_ != oneSidedDef || bestFitR < initialR) { 
        // must do constrained fit (if there's something to fit besides XS)
        //std::cout << "PERFORMING CONSTRAINED FIT " << r->GetName() << " == " << r->getVal() << std::endl;
        thisNLL = (nuisances_.getSize() > 0 ? minNLL(/*constrained=*/true, r) : nll_->getVal());
        if (thisNLL - nullNLL < -0.02) { 
            DBG(DBG_PLTestStat_main, (printf("  --> constrained fit is better... will repeat unconstrained fit\n")))
            utils::setAllConstant(poiParams_,false);
            nullNLL = minNLL(/*constrained=*/false, r);
            bestFitR = r->getVal();
            if (bestFitR > initialR && oneSided_ == oneSidedDef) {
                DBG(DBG_PLTestStat_main, (printf("   after re-fit, signal %7.4f > %7.4f, test statistics will be zero.\n", bestFitR, initialR)))
                thisNLL = nullNLL;
            }
        }
        /* This piece of debug code was added to see if we were finding a local minimum at zero instead of the global minimum.
           It doesn't seem to be the case, however  
        if (do_rescan && fabs(thisNLL - nullNLL) < 0.2 && initialR == 0) {
            if (do_debug) printf("Doing full rescan. best fit r = %.4f, -lnQ = %.5f\n", bestFitR, thisNLL-nullNLL);
            for (double rx = 0; rx <= 5; rx += 0.01) {
                r->setVal(rx); double y = nll_->getVal();
                if (y < nullNLL) {
                    printf("%4.2f\t%8.5f\t<<==== ALERT\n", rx, y - nullNLL);
                } else {
                    printf("%4.2f\t%8.5f\n", rx, y - nullNLL);
                }
            }
        } */
        if (bestFitR > initialR && oneSided_ == signFlipDef) {
            DBG(DBG_PLTestStat_main, (printf("   fitted signal %7.4f > %7.4f, test statistics will be negative.\n", bestFitR, initialR)))
            std::swap(thisNLL, nullNLL);
        }
    } else {
        DBG(DBG_PLTestStat_main, (printf("   signal fit %7.4f > %7.4f: don't need to compute numerator\n", bestFitR, initialR)))
    }

    //Restore initial state, to avoid issues with ToyMCSampler
    *params_ = initialState;

    if (!canKeepNLL) nll_.reset();

    DBG(DBG_PLTestStat_main, (printf("\nNLLs:  num % 10.4f, den % 10.4f (signal %7.4f), test stat % 10.4f\n", thisNLL, nullNLL, bestFitR, thisNLL-nullNLL)))
    if (do_debug) printf("\nNLLs:  num % 10.4f, den % 10.4f (signal %7.4f), test stat % 10.4f\n", thisNLL, nullNLL, bestFitR, thisNLL-nullNLL);
    //return std::min(thisNLL-nullNLL, 0.);
    return thisNLL-nullNLL;
}

virtual const TString ProfiledLikelihoodTestStatOpt::GetVarName ( ) const [inline, virtual]

Definition at line 61 of file ProfiledLikelihoodRatioTestStatExt.h.

{ return "- log (#lambda)"; }

double ProfiledLikelihoodTestStatOpt::minNLL	(	bool	constrained,
		RooRealVar *	r = `0`
	)		`[private]`

Definition at line 374 of file ProfiledLikelihoodRatioTestStatExt.cc.

References CascadeMinimizer::Constrained, CascadeMinimizer::minimize(), mode, nll_, CascadeMinimizer::Unconstrained, and verbosity_.

Referenced by Evaluate().

{
    CascadeMinimizer::Mode mode(constrained ? CascadeMinimizer::Constrained : CascadeMinimizer::Unconstrained);
    CascadeMinimizer minim(*nll_, mode, r);
    minim.minimize(verbosity_-2);
    return nll_->getVal();
}