Save L-curve and other curves, as well as solutions.
352 cout << __LINE__ <<
"\tScanning L-curve" << endl;
354 TGraph * lCurve =
nullptr;
355 TSpline * logTauX =
nullptr,
357 * logTauCurvature =
nullptr;
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);
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);
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;
381 double tau = density->GetTau();
382 cout << __LINE__ <<
"\ttau = " << tau << endl;
384 lCurve ->SetNameTitle(
"lCurve" , Form(
"#tau = %f", tau));
385 bestLcurve->SetNameTitle(
"bestLcurve", Form(
"%f", tau));
387 logTauX ->SetNameTitle(
"logTauX",
"log(#chi_{A}^{2})");
388 bestLogTauX->SetNameTitle(
"bestLogTauX", Form(
"#chi_{A}^{2}/ndf = %f", chi2a/ndf));
390 logTauY ->SetNameTitle(
"logTauY",
"log(#chi_{L}^{2}/#tau)");
391 bestLogTauY->SetNameTitle(
"bestLogTauY", Form(
"#chi_{L}^{2}/ndf = %f", chi2l/ndf));
393 logTauCurvature ->SetNameTitle(
"logTauCurvature", Form(
"%f", tau));
394 bestLogTauCurvature->SetNameTitle(
"bestLogTauCurvature", Form(
"%f", tau));
397 lCurve->Write(); bestLcurve->Write();
398 logTauX->Write(); bestLogTauX->Write();
399 logTauY->Write(); bestLogTauY->Write();
400 logTauCurvature->Write(); bestLogTauCurvature->Write();
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);
412 delete logTauCurvature;