DAS::Unfolding::Tikhonov Namespace Reference

Functions
void	PrepareRegularisation (TString MCvar, const char *fname, unique_ptr< MyTUnfoldDensity > &density)
void	RedoUnfold (unique_ptr< MyTUnfoldDensity > &density, double newtau, TString name, TString subtitle, shared_ptr< TFile > fOut, const map< TString, unique_ptr< TH1 >> &miss)
double	ScanLcurve (unique_ptr< MyTUnfoldDensity > &density, TDirectory *dOut, const pt::ptree &config)

Function Documentation

PrepareRegularisation()

void DAS::Unfolding::Tikhonov::PrepareRegularisation	(	TString	MCvar,
		const char *	fname,
		unique_ptr< MyTUnfoldDensity > &	density
	)

Prepare regularisation

Load bias and L-matrix from the output file of one of the getLmatrix* commands

Todo:

calculate \( A^T_{yy} V^{-1} A \) & adapt the L matrix to have the same rank

{
    cout << __LINE__ << '\t' << __func__ << endl;
    cerr << __LINE__ << '\t' << "Warning: the implementation of the regularisation "
                                "is not yet stable. For instance, it is not expected "
                                "to work correctly in case of empty bins. Use at your "
                                "own risk.\n" << def;
 
    auto fReg = make_unique<TFile>(fname, "READ");
 
    cout << __LINE__ << "\tgetting bias" << endl;
    auto bias = unique_ptr<TH1>(fReg->Get<TH1>(MCvar + "/bias"));
    if (!bias) BOOST_THROW_EXCEPTION( DE::BadInput("Could not find the bias", fReg) );
    density->SetBias(bias.get());
 
    // getting L matrix
    cout << __LINE__ << "\tgetting L-matrix" << endl;
    auto L = unique_ptr<TH2>(fReg->Get<TH2>(MCvar + "/L"));
    if (!L) BOOST_THROW_EXCEPTION( DE::BadInput("Could not find the L-matrix", fReg) );
 
    int Nconditions = L->GetNbinsY(), // arbitrarily large number
        NtrueBins   = L->GetNbinsX(); // must be exactly the number of bins at particle level
 
 
    cout << __LINE__ << "\tlooping over " << Nconditions << " conditions" << endl;
    for (int j = 1; j <= Nconditions; ++j) { // y-axis
 
        vector<int> indices;
        vector<double> conditions;
        indices.reserve(Nconditions);
        conditions.reserve(Nconditions);
 
        for (int i = 1; i <= NtrueBins; ++i) { // x-axis
            double content = L->GetBinContent(i, j);
            if (std::abs(content) < deps) continue;
            indices.push_back(i);
            conditions.push_back(content);
        }
 
        // cheating with call to protected method
        if (indices.size() != conditions.size())
            BOOST_THROW_EXCEPTION( DE::BadInput("inconsistent number of indices and conditions",
                                                fReg) );
 
        density->MyAddRegularisationCondition(indices.size(), indices.data(), conditions.data());
    }
    cout << __LINE__ << "\tregularisation is ready" << endl;
}

RedoUnfold()

void DAS::Unfolding::Tikhonov::RedoUnfold	(	unique_ptr< MyTUnfoldDensity > &	density,
		double	newtau,
		TString	name,
		TString	subtitle,
		shared_ptr< TFile >	fOut,
		const map< TString, unique_ptr< TH1 >> &	miss
	)

RedoUnfold

Repeat unfolding for alternative values of tau parameter

Todo:

cov (needed for smoothing)

{
    auto dOut = fOut->mkdir(name, subtitle);
    density->DoUnfold(newtau);
    auto unf = unique_ptr<TH1>(density->GetOutput("unfold", "particle level (" + subtitle + ")"));
    CorrectInefficiencies(unf, miss);
    Checks::NegativeBins(unf);
    unf->SetDirectory(dOut);
    dOut->cd();
    unf->GetXaxis()->SetTitle("gen");
    unf->Write();
}

ScanLcurve()

double DAS::Unfolding::Tikhonov::ScanLcurve	(	unique_ptr< MyTUnfoldDensity > &	density,
		TDirectory *	dOut,
		const pt::ptree &	config
	)

Save L-curve and other curves, as well as solutions.

{
    cout << __LINE__ << "\tScanning L-curve" << endl;
 
    TGraph * lCurve = nullptr;
    TSpline * logTauX = nullptr,
            * logTauY = nullptr,
            * logTauCurvature = nullptr;
 
    auto nPoints = config.get_optional<int>  ("unfolding.TUnfold.regularisation.Lcurve.nPoints").value_or(100);
    auto mintau = config.get_optional<double>("unfolding.TUnfold.regularisation.Lcurve.mintau").value_or(1e-6),
         maxtau = config.get_optional<double>("unfolding.TUnfold.regularisation.Lcurve.maxtau").value_or(1e+2);
    int iBest = density->ScanLcurve(nPoints, mintau, maxtau,
                                    &lCurve, &logTauX, &logTauY, &logTauCurvature);
 
    cout << __LINE__ << "\tSaving control plots for L-curve scan" << endl;
    Double_t t[1],x[1],y[1], c[1];
    logTauX->GetKnot(iBest,t[0],x[0]);
    logTauY->GetKnot(iBest,t[0],y[0]);
    logTauCurvature->GetKnot(iBest,t[0],c[0]);
    auto bestLcurve          = make_unique<TGraph>(1,x,y);
    auto bestLogTauX         = make_unique<TGraph>(1,t,x);
    auto bestLogTauY         = make_unique<TGraph>(1,t,y);
    auto bestLogTauCurvature = make_unique<TGraph>(1,t,c);
 
    double chi2a = density->GetChi2A(),
           chi2l = density->GetChi2L();
    int ndf = density->GetNdf();
    if (ndf == 0) BOOST_THROW_EXCEPTION( runtime_error("Can't regularise if ndf = 0") );
    cout << __LINE__ << "\tchi^2/ndf = " << chi2a/ndf << " + " << chi2l/ndf << endl;
 
    double tau = density->GetTau();
    cout << __LINE__ << "\ttau = " << tau << endl;
 
    lCurve    ->SetNameTitle("lCurve"    , Form("#tau = %f", tau));
    bestLcurve->SetNameTitle("bestLcurve", Form(       "%f", tau));
 
    logTauX    ->SetNameTitle("logTauX", "log(#chi_{A}^{2})");
    bestLogTauX->SetNameTitle("bestLogTauX", Form("#chi_{A}^{2}/ndf = %f", chi2a/ndf));
 
    logTauY    ->SetNameTitle("logTauY", "log(#chi_{L}^{2}/#tau)");
    bestLogTauY->SetNameTitle("bestLogTauY", Form("#chi_{L}^{2}/ndf = %f", chi2l/ndf));
 
    logTauCurvature    ->SetNameTitle("logTauCurvature", Form("%f", tau));
    bestLogTauCurvature->SetNameTitle("bestLogTauCurvature", Form("%f", tau));
 
    dOut->cd();
    lCurve->Write();           bestLcurve->Write();
    logTauX->Write();          bestLogTauX->Write();
    logTauY->Write();          bestLogTauY->Write();
    logTauCurvature->Write();  bestLogTauCurvature->Write();
 
    auto L  = unique_ptr<TH2>(density->GetL          ("L"          , "L-matrix"         ));
    auto Lx = unique_ptr<TH1>(density->GetLxMinusBias("LxMinusBias", "L*(x-f_{b}*x_{0})"));
    L ->SetDirectory(dOut);
    Lx->SetDirectory(dOut);
    L ->Write();
    Lx->Write();
 
    delete lCurve;
    delete logTauX;
    delete logTauY;
    delete logTauCurvature;
 
    return tau;
}