DAS  3.0
Das Analysis System
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, unique_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
277 {
278  cout << __LINE__ << '\t' << __func__ << endl;
279  cerr << __LINE__ << '\t' << "Warning: the implementation of the regularisation "
280  "is not yet stable. For instance, it is not expected "
281  "to work correctly in case of empty bins. Use at your "
282  "own risk.\n" << def;
283 
284  auto fReg = make_unique<TFile>(fname, "READ");
285 
286  cout << __LINE__ << "\tgetting bias" << endl;
287  auto bias = unique_ptr<TH1>(fReg->Get<TH1>(MCvar + "/bias"));
288  if (!bias) BOOST_THROW_EXCEPTION( DE::BadInput("Could not find the bias", fReg) );
289  density->SetBias(bias.get());
290 
291  // getting L matrix
292  cout << __LINE__ << "\tgetting L-matrix" << endl;
293  auto L = unique_ptr<TH2>(fReg->Get<TH2>(MCvar + "/L"));
294  if (!L) BOOST_THROW_EXCEPTION( DE::BadInput("Could not find the L-matrix", fReg) );
295 
296  int Nconditions = L->GetNbinsY(), // arbitrarily large number
297  NtrueBins = L->GetNbinsX(); // must be exactly the number of bins at particle level
298 
300 
301  cout << __LINE__ << "\tlooping over " << Nconditions << " conditions" << endl;
302  for (int j = 1; j <= Nconditions; ++j) { // y-axis
303 
304  vector<int> indices;
305  vector<double> conditions;
306  indices.reserve(Nconditions);
307  conditions.reserve(Nconditions);
308 
309  for (int i = 1; i <= NtrueBins; ++i) { // x-axis
310  double content = L->GetBinContent(i, j);
311  if (std::abs(content) < deps) continue;
312  indices.push_back(i);
313  conditions.push_back(content);
314  }
315 
316  // cheating with call to protected method
317  if (indices.size() != conditions.size())
318  BOOST_THROW_EXCEPTION( DE::BadInput("inconsistent number of indices and conditions",
319  fReg) );
320 
321  density->MyAddRegularisationCondition(indices.size(), indices.data(), conditions.data());
322  }
323  cout << __LINE__ << "\tregularisation is ready" << endl;
324 }

◆ RedoUnfold()

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

RedoUnfold

Repeat unfolding for alternative values of tau parameter

Todo:
cov (needed for smoothing)
333 {
334  auto dOut = fOut->mkdir(name, subtitle);
335  density->DoUnfold(newtau);
336  auto unf = unique_ptr<TH1>(density->GetOutput("unfold", "particle level (" + subtitle + ")"));
337  CorrectInefficiencies(unf, miss);
339  unf->SetDirectory(dOut);
340  dOut->cd();
341  unf->GetXaxis()->SetTitle("gen");
342  unf->Write();
344 }

◆ 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.

351 {
352  cout << __LINE__ << "\tScanning L-curve" << endl;
353 
354  TGraph * lCurve = nullptr;
355  TSpline * logTauX = nullptr,
356  * logTauY = nullptr,
357  * logTauCurvature = nullptr;
358 
359  auto nPoints = config.get_optional<int> ("unfolding.TUnfold.regularisation.Lcurve.nPoints").value_or(100);
360  auto mintau = config.get_optional<double>("unfolding.TUnfold.regularisation.Lcurve.mintau").value_or(1e-6),
361  maxtau = config.get_optional<double>("unfolding.TUnfold.regularisation.Lcurve.maxtau").value_or(1e+2);
362  int iBest = density->ScanLcurve(nPoints, mintau, maxtau,
363  &lCurve, &logTauX, &logTauY, &logTauCurvature);
364 
365  cout << __LINE__ << "\tSaving control plots for L-curve scan" << endl;
366  Double_t t[1],x[1],y[1], c[1];
367  logTauX->GetKnot(iBest,t[0],x[0]);
368  logTauY->GetKnot(iBest,t[0],y[0]);
369  logTauCurvature->GetKnot(iBest,t[0],c[0]);
370  auto bestLcurve = make_unique<TGraph>(1,x,y);
371  auto bestLogTauX = make_unique<TGraph>(1,t,x);
372  auto bestLogTauY = make_unique<TGraph>(1,t,y);
373  auto bestLogTauCurvature = make_unique<TGraph>(1,t,c);
374 
375  double chi2a = density->GetChi2A(),
376  chi2l = density->GetChi2L();
377  int ndf = density->GetNdf();
378  if (ndf == 0) BOOST_THROW_EXCEPTION( runtime_error("Can't regularise if ndf = 0") );
379  cout << __LINE__ << "\tchi^2/ndf = " << chi2a/ndf << " + " << chi2l/ndf << endl;
380 
381  double tau = density->GetTau();
382  cout << __LINE__ << "\ttau = " << tau << endl;
383 
384  lCurve ->SetNameTitle("lCurve" , Form("#tau = %f", tau));
385  bestLcurve->SetNameTitle("bestLcurve", Form( "%f", tau));
386 
387  logTauX ->SetNameTitle("logTauX", "log(#chi_{A}^{2})");
388  bestLogTauX->SetNameTitle("bestLogTauX", Form("#chi_{A}^{2}/ndf = %f", chi2a/ndf));
389 
390  logTauY ->SetNameTitle("logTauY", "log(#chi_{L}^{2}/#tau)");
391  bestLogTauY->SetNameTitle("bestLogTauY", Form("#chi_{L}^{2}/ndf = %f", chi2l/ndf));
392 
393  logTauCurvature ->SetNameTitle("logTauCurvature", Form("%f", tau));
394  bestLogTauCurvature->SetNameTitle("bestLogTauCurvature", Form("%f", tau));
395 
396  dOut->cd();
397  lCurve->Write(); bestLcurve->Write();
398  logTauX->Write(); bestLogTauX->Write();
399  logTauY->Write(); bestLogTauY->Write();
400  logTauCurvature->Write(); bestLogTauCurvature->Write();
401 
402  auto L = unique_ptr<TH2>(density->GetL ("L" , "L-matrix" ));
403  auto Lx = unique_ptr<TH1>(density->GetLxMinusBias("LxMinusBias", "L*(x-f_{b}*x_{0})"));
404  L ->SetDirectory(dOut);
405  Lx->SetDirectory(dOut);
406  L ->Write();
407  Lx->Write();
408 
409  delete lCurve;
410  delete logTauX;
411  delete logTauY;
412  delete logTauCurvature;
413 
414  return tau;
415 }
DYToLL_M-50_13TeV_pythia8_cff_GEN_SIM_RECOBEFMIX_DIGI_L1_DIGI2RAW_L1Reco_RECO.name
name
Definition: DYToLL_M-50_13TeV_pythia8_cff_GEN_SIM_RECOBEFMIX_DIGI_L1_DIGI2RAW_L1Reco_RECO.py:48
Ntupliser_cfg.cerr
cerr
Definition: Ntupliser_cfg.py:93
jercExample.fname
fname
Definition: jercExample.py:88
DAS::Unfolding::Checks::NegativeBins
void NegativeBins(const unique_ptr< TH1 > &unf)
Definition: unfold.cc:122
Step::def
static const char * def
Definition: Step.h:36
deps
static const auto deps
Definition: unfold.cc:28
DAS::Unfolding::CorrectInefficiencies
void CorrectInefficiencies(unique_ptr< TH1 > &unf, const map< TString, unique_ptr< TH1 >> &miss)
Definition: unfold.cc:58
Ntupliser_cfg.config
config
Definition: Ntupliser_cfg.py:263
Darwin::Exceptions::BadInput
Generic exception for ill-defined input (before the event loop).
Definition: exceptions.h:68