DAS::Unfolding::Checks Namespace Reference

Functions
void	NegativeBins (const unique_ptr< TH1 > &unf)
unique_ptr< TH2 >	GetPM (const unique_ptr< TH2 > &RM)
void	folding (const unique_ptr< TH1 > &genIn, const unique_ptr< TH1 > &recIn, const unique_ptr< TH2 > &PMin, const map< TString, unique_ptr< TH1 >> &miss, const map< TString, unique_ptr< TH1 >> &fake)
void	BLT (const unique_ptr< TH1 > &dataDist, const unique_ptr< TH2 > &dataCov, const unique_ptr< TH1 > &MC, int rebin=1)
void	CovMargin (unique_ptr< TH1 > &h, unique_ptr< TH2 > &cov, TString name, int ibin)

Function Documentation

BLT()

void DAS::Unfolding::Checks::BLT	(	const unique_ptr< TH1 > &	dataDist,
		const unique_ptr< TH2 > &	dataCov,
		const unique_ptr< TH1 > &	MC,
		int	rebin = 1
	)

Bottom Line Test

Checks the agreement of data and MC before and after unfolding

Todo:

current implementation is insufficient for complex binning schemes

{
    cout << __LINE__ << '\t' << __func__ << endl;
    cerr << __LINE__ << '\t' << orange << "Warning: the BLT only gives reasonable "
                                          "results for simples rebinning schemes "
                                          "(e.g. 2 rec bins for 1 gen bins)\n" << def;
    auto reb_dataDist = Clone<TH1>(dataDist, "reb_dataDist");
    auto reb_dataCov  = Clone<TH2>(dataCov,  "reb_dataCov" );
    auto reb_MC       = Clone<TH1>(MC,       "reb_MC"      );
 
    reb_dataDist->Rebin(rebin);
    reb_dataCov->Rebin2D(rebin, rebin);
    reb_MC->Rebin(rebin);
    vector<int> indices;
    for (int i = 1; i <= dataCov->GetNbinsX(); ++i)
        if (dataCov->GetBinContent(i,i) > 0) indices.push_back(i);
    int ndf = indices.size();
 
    TVectorD v(ndf);
    TMatrixD M(ndf,ndf);
    M.SetTol(1e-100);
    for (int i = 0; i < ndf; ++i) {
        int index = indices.at(i);
        double iData = dataDist->GetBinContent(index),
               iMC   = MC      ->GetBinContent(index);
        v(i) = iData-iMC;
        for (int j = 0; j < ndf; ++j) {
            int jndex = indices.at(j);
            M(i,j) = dataCov->GetBinContent(index,jndex);
        }
    }
    M.Invert();
 
    double chi2 = v*(M*v);
 
    cout << "chi2/ndf = " << chi2 << " / " << ndf << " = " << (chi2/ndf) << endl;
}

CovMargin()

void DAS::Unfolding::Checks::CovMargin	(	unique_ptr< TH1 > &	h,
		unique_ptr< TH2 > &	cov,
		TString	name,
		int	ibin
	)

Set underflow/overflow to 0.

{
    if (   h->GetBinContent(ibin) == 0
        && h->GetBinError(ibin) == 0
        && cov->GetBinContent(ibin,ibin) == 0) return;
 
    cerr << __LINE__ << orange << "\tWarning: " << name << " is not empty! "
            "It will be set to zero, but no guarantee it works...\n" << def;
 
    h->SetBinContent(ibin,0);
    h->SetBinError  (ibin,0);
    for (int i = 0; i <= h->GetNbinsX(); ++i) {
        cov->SetBinContent(ibin,i,0);
        cov->SetBinContent(i,ibin,0);
    }
}

folding()

void DAS::Unfolding::Checks::folding	(	const unique_ptr< TH1 > &	genIn,
		const unique_ptr< TH1 > &	recIn,
		const unique_ptr< TH2 > &	PMin,
		const map< TString, unique_ptr< TH1 >> &	miss,
		const map< TString, unique_ptr< TH1 >> &	fake
	)

Folding test

Check that the MC closes: (gen - sum of misses) * PM = (rec - sum of fakes)

Parameters

genIn	input gen-level spectrum
recIn	input rec-level spectrum
PMin	input probability matrix
miss	miss counts
fake	fake counts

{
    cout << __LINE__ << '\t' << __func__ << endl;
 
    unique_ptr<TH1> gen = Clone<TH1>(genIn, "gen_folding");
    unique_ptr<TH1> rec = Clone<TH1>(recIn, "rec_folding");
    unique_ptr<TH2> PM  = Clone<TH2>(PMin ,  "PM_folding");
 
    for (auto& m: miss) gen->Add(m.second.get(), -1);
    for (auto& f: fake) rec->Add(f.second.get(), -1);
 
    const int Nx = PM->GetNbinsX(),
              Ny = PM->GetNbinsY();
 
    // check folding
    bool Fail = false;
    for (int j = 1; j <= Ny; ++j) {
        double error = rec->GetBinError(j);
        if (::abs(error) < deps) continue;
 
        double content = rec->GetBinContent(j);
        Fail |= isnan(content);
 
        for (int i = 1; i <= Nx; ++i)
            content -= gen->GetBinContent(i) * PM->GetBinContent(i,j);
 
        if (::abs(content) < error / 100) continue;
        Fail = true;
        cerr << __LINE__ << orange << "\tWarning: Folding test is failing in bin "
             << j << " with (gen - miss) x PM - (rec - fake) = " << content
             << " > error / 100 = " << (error / 100) << def << '\n';
    }
 
    if (Fail) /*BOOST_THROW_EXCEPTION( std::runtime_error(*/ cerr << red << "Failure: folding test has failed\n" << def /*) )*/;
}

GetPM()

unique_ptr<TH2> DAS::Unfolding::Checks::GetPM

(

const unique_ptr< TH2 > &

RM

)

Calculate probability matrix.

{
    unique_ptr<TH2> PM = Clone<TH2>(RM, "PM");
 
    const int Nx = RM->GetNbinsX(),
              Ny = RM->GetNbinsY();
 
    // normalisation
    for (int i = 1; i <= Nx; ++i) {
        double normalisation = RM->Integral(i, i, 1, Ny);
        if (normalisation > 0)
        for (int j = 1; j <= Ny; ++j) {
            double content = RM->GetBinContent(i,j);
            content /= normalisation;
            PM->SetBinContent(i,j, content);
        }
    }
 
    return PM;
}

NegativeBins()

void DAS::Unfolding::Checks::NegativeBins

(

const unique_ptr< TH1 > &

unf

)

Check negative entries

Note: you should not necessarily worry if you get a value compatible with zero BUT you should not trust the result if the value is not at all compatible with zero

{
    cout << __LINE__ << '\t' << __func__ << endl;
    for (int i = 1; i <= unf->GetNbinsX(); ++i) {
        double content = unf->GetBinContent(i);
        if (content >= 0) continue;
        double error = unf->GetBinError(i);
        cerr << __LINE__ << orange << "\tWarning: bin " << i
             << " has negative content: " << content << " +/- " << error << '\n' << def;
    }
}