Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

refactor: residuals and pulls analysis #1892

Merged
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 5 additions & 2 deletions Examples/Scripts/TrackingPerformance/ResidualsAndPulls.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,9 @@ int main(int argc, char** argv) {
ao("predicted", bool_switch(), "Analyze the predicted parameters.");
ao("filtered", bool_switch(), "Analyze the filtered parameters.");
ao("smoothed", bool_switch(), "Analyze the smoothed parameters.");
ao("fit", bool_switch(), "Fit the smoothed parameters.");
ao("fit-predicted", bool_switch(), "Fit the predicted parameters.");
ao("fit-filtered", bool_switch(), "Fit the filtered parameters.");
ao("fit-smoothed", bool_switch(), "Fit the smoothed parameters.");
ao("save", value<std::string>()->default_value("png"),
"Output save format (to be interpreted by ROOT).");

Expand Down Expand Up @@ -66,7 +68,8 @@ int main(int argc, char** argv) {
switch (boundParamResolution(
iFile, iTree, oFile, vm["predicted"].as<bool>(),
vm["filtered"].as<bool>(), vm["smoothed"].as<bool>(),
vm["fit"].as<bool>(), saveAs)) {
vm["fit-predicted"].as<bool>(), vm["fit-filtered"].as<bool>(),
vm["fit-smoothed"].as<bool>(), saveAs)) {
case -1: {
std::cout << "*** Input file could not be opened, check name/path."
<< std::endl;
Expand Down
175 changes: 99 additions & 76 deletions Examples/Scripts/TrackingPerformance/boundParamResolution.C
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,7 @@ using namespace ROOT;
int boundParamResolution(const std::string& inFile, const std::string& treeName,
const std::string& outFile, bool predicted = true,
bool filtered = true, bool smoothed = true,
bool fitFiltered = true, bool fitPredicted = true,
bool fitSmoothed = true,
const std::string& saveAs = "") {
// Some style options
Expand Down Expand Up @@ -411,10 +412,10 @@ int boundParamResolution(const std::string& inFile, const std::string& treeName,
for (unsigned int i = 0; i < tsReader.nMeasurements; i++) {
// global profile filling
if (predicted and tsReader.predicted->at(i)) {
float x_prt = tsReader.g_x_prt->at(i);
float y_prt = tsReader.g_y_prt->at(i);
float r_prt = std::sqrt(x_prt * x_prt + y_prt * y_prt);
float z_prt = tsReader.g_z_prt->at(i);
auto x_prt = tsReader.g_x_prt->at(i);
auto y_prt = tsReader.g_y_prt->at(i);
auto r_prt = std::sqrt(x_prt * x_prt + y_prt * y_prt);
auto z_prt = tsReader.g_z_prt->at(i);
p2d_res_zr_prt[0]->Fill(z_prt, r_prt, tsReader.res_LOC0_prt->at(i));
p2d_res_zr_prt[1]->Fill(z_prt, r_prt, tsReader.res_LOC1_prt->at(i));
p2d_res_zr_prt[2]->Fill(z_prt, r_prt, tsReader.res_PHI_prt->at(i));
Expand All @@ -429,10 +430,10 @@ int boundParamResolution(const std::string& inFile, const std::string& treeName,
p2d_pull_zr_prt[5]->Fill(z_prt, r_prt, tsReader.pull_T_prt->at(i));
}
if (filtered and tsReader.filtered->at(i)) {
float x_flt = tsReader.g_x_flt->at(i);
float y_flt = tsReader.g_y_flt->at(i);
float r_flt = std::sqrt(x_flt * x_flt + y_flt * y_flt);
float z_flt = tsReader.g_z_flt->at(i);
auto x_flt = tsReader.g_x_flt->at(i);
auto y_flt = tsReader.g_y_flt->at(i);
auto r_flt = std::sqrt(x_flt * x_flt + y_flt * y_flt);
auto z_flt = tsReader.g_z_flt->at(i);
p2d_res_zr_flt[0]->Fill(z_flt, r_flt, tsReader.res_LOC0_flt->at(i));
p2d_res_zr_flt[1]->Fill(z_flt, r_flt, tsReader.res_LOC1_flt->at(i));
p2d_res_zr_flt[2]->Fill(z_flt, r_flt, tsReader.res_PHI_flt->at(i));
Expand All @@ -447,10 +448,10 @@ int boundParamResolution(const std::string& inFile, const std::string& treeName,
p2d_pull_zr_flt[5]->Fill(z_flt, r_flt, tsReader.pull_T_flt->at(i));
}
if (smoothed and tsReader.smoothed->at(i)) {
float x_smt = tsReader.g_x_smt->at(i);
float y_smt = tsReader.g_y_smt->at(i);
float r_smt = std::sqrt(x_smt * x_smt + y_smt * y_smt);
float z_smt = tsReader.g_z_smt->at(i);
auto x_smt = tsReader.g_x_smt->at(i);
auto y_smt = tsReader.g_y_smt->at(i);
auto r_smt = std::sqrt(x_smt * x_smt + y_smt * y_smt);
auto z_smt = tsReader.g_z_smt->at(i);
p2d_res_zr_smt[0]->Fill(z_smt, r_smt, tsReader.res_LOC0_smt->at(i));
p2d_res_zr_smt[1]->Fill(z_smt, r_smt, tsReader.res_LOC1_smt->at(i));
p2d_res_zr_smt[2]->Fill(z_smt, r_smt, tsReader.res_PHI_smt->at(i));
Expand Down Expand Up @@ -535,12 +536,12 @@ int boundParamResolution(const std::string& inFile, const std::string& treeName,
// Section 5: Histogram plotting

// Plotting global profiles
TCanvas* respull_mean_prf =
auto respull_mean_prf =
new TCanvas("respull_mean_prf",
"Residual/Pull Distributions - mean profiles", 1800, 800);
respull_mean_prf->Divide(3, 2);

TCanvas* respull_var_prf =
auto respull_var_prf =
new TCanvas("respull_var_prf",
"Residual/Pull Distributions - variance profiles", 1800, 800);
respull_var_prf->Divide(3, 2);
Expand Down Expand Up @@ -573,8 +574,8 @@ int boundParamResolution(const std::string& inFile, const std::string& treeName,
respull_var_prf->cd(ipar + 1);
auto zAxis = profiles[ipar]->GetXaxis();
auto rAxis = profiles[ipar]->GetYaxis();
int binsZ = zAxis->GetNbins();
int binsR = rAxis->GetNbins();
auto binsZ = zAxis->GetNbins();
auto binsR = rAxis->GetNbins();
std::string hist_name = "all_";
hist_name += res_pull;
hist_name += "_";
Expand Down Expand Up @@ -624,11 +625,11 @@ int boundParamResolution(const std::string& inFile, const std::string& treeName,
}

// Plotting residual
TCanvas* residuals =
auto residuals =
new TCanvas("residuals", "Residual Distributions", 1200, 800);
residuals->Divide(3, 2);

TCanvas* pulls = new TCanvas("pulls", "Pull distributions", 1200, 800);
auto pulls = new TCanvas("pulls", "Pull distributions", 1200, 800);
pulls->Divide(3, 2);

for (auto [vol, layers] : volLayIds) {
Expand All @@ -639,30 +640,27 @@ int boundParamResolution(const std::string& inFile, const std::string& treeName,
for (size_t ipar = 0; ipar < paramNames.size(); ipar++) {
residuals->cd(ipar + 1);

TLegend* legend = new TLegend(0.7, 0.7, 0.9, 0.9);
auto legend = new TLegend(0.7, 0.7, 0.9, 0.9);

const std::string name = vlID + paramNames.at(ipar);

// Draw them
if (smoothed) {
res_smt[vlID + paramNames.at(ipar)]->DrawNormalized("");
res_smt[vlID + paramNames.at(ipar)]->Write();
legend->AddEntry(res_smt[vlID + paramNames.at(ipar)], "smoothed",
"l");
res_smt[name]->DrawNormalized("");
res_smt[name]->Write();
legend->AddEntry(res_smt[name], "smoothed", "l");
}
if (filtered) {
std::string drawOptions = smoothed ? "same" : "";
res_flt[vlID + paramNames.at(ipar)]->DrawNormalized(
drawOptions.c_str());
res_flt[vlID + paramNames.at(ipar)]->Write();
legend->AddEntry(res_flt[vlID + paramNames.at(ipar)], "filtered",
"l");
res_flt[name]->DrawNormalized(drawOptions.c_str());
res_flt[name]->Write();
legend->AddEntry(res_flt[name], "filtered", "l");
}
if (predicted) {
std::string drawOptions = (smoothed or filtered) ? "same" : "";
res_prt[vlID + paramNames.at(ipar)]->DrawNormalized(
drawOptions.c_str());
res_prt[vlID + paramNames.at(ipar)]->Write();
legend->AddEntry(res_prt[vlID + paramNames.at(ipar)], "predicted",
"l");
res_prt[name]->DrawNormalized(drawOptions.c_str());
res_prt[name]->Write();
legend->AddEntry(res_prt[name], "predicted", "l");
}

legend->SetBorderSize(0);
Expand All @@ -676,69 +674,94 @@ int boundParamResolution(const std::string& inFile, const std::string& treeName,

// Pull plotting & writing
for (size_t ipar = 0; ipar < paramNames.size(); ipar++) {
const std::string name = vlID + paramNames.at(ipar);

pulls->cd(ipar + 1);

TLegend* legend = new TLegend(0.7, 0.7, 0.9, 0.9);
auto legend = new TLegend(0.7, 0.5, 0.9, 0.9);

// Fit the smoothed distribution
if (smoothed) {
Double_t scale =
1. / pull_smt[vlID + paramNames.at(ipar)]->Integral();
pull_smt[vlID + paramNames.at(ipar)]->Sumw2();
pull_smt[vlID + paramNames.at(ipar)]->Scale(scale);
pull_smt[vlID + paramNames.at(ipar)]->Draw("pe");
pull_smt[vlID + paramNames.at(ipar)]->Write();
auto drawOptions = "pe";

legend->AddEntry(pull_smt[vlID + paramNames.at(ipar)], "smoothed",
"pe");
auto scale = 1. / pull_smt[name]->Integral("width");
pull_smt[name]->Scale(scale);
pull_smt[name]->Draw(drawOptions);
pull_smt[name]->Write();

legend->AddEntry(pull_smt[name], "smoothed", "pe");

if (fitSmoothed) {
pull_smt[vlID + paramNames.at(ipar)]->Fit("gaus", "q");
TF1* gauss =
pull_smt[vlID + paramNames.at(ipar)]->GetFunction("gaus");
gauss->SetLineColor(kGreen);
float mu = gauss->GetParameter(1);
float sigma = gauss->GetParameter(2);

// Draw the sigma
TString mu_info;
mu_info += mu;
TString sigma_info;
sigma_info += sigma;

TString mfit_info = "#mu = ";
mfit_info += mu_info(0, 5);

TString sfit_info = "#sigma = ";
sfit_info += sigma_info(0, 5);

legend->AddEntry(gauss, sfit_info.Data(), "l");
legend->AddEntry(pull_smt[vlID + paramNames.at(ipar)],
mfit_info.Data(), "");
pull_smt[name]->Fit("gaus", "q");
TF1* gauss = pull_smt[name]->GetFunction("gaus");
gauss->SetLineColorAlpha(kBlack, 0.5);

auto mu = gauss->GetParameter(1);
auto sigma = gauss->GetParameter(2);
auto mu_fit_info = TString::Format("#mu = %.3f", mu);
auto su_fit_info = TString::Format("#sigma = %.3f", sigma);

legend->AddEntry(gauss, mu_fit_info.Data(), "l");
legend->AddEntry(pull_prt[name], su_fit_info.Data(), "");
}
}

if (filtered) {
std::string drawOptions = smoothed ? "same" : "";
auto drawOptions = smoothed ? "pe same" : "pe";

auto scale = 1. / pull_flt[name]->Integral("width");
pull_flt[name]->Scale(scale);
pull_flt[name]->Draw(drawOptions);
pull_flt[name]->Write();

legend->AddEntry(pull_flt[name], "filtered", "pe");

pull_flt[vlID + paramNames.at(ipar)]->DrawNormalized(
drawOptions.c_str());
pull_flt[vlID + paramNames.at(ipar)]->Write();
if (fitFiltered) {
pull_flt[name]->Fit("gaus", "q", "same");
TF1* gauss = pull_flt[name]->GetFunction("gaus");
gauss->SetLineColorAlpha(kBlue, 0.5);

legend->AddEntry(pull_flt[vlID + paramNames.at(ipar)], "filtered",
"lp");
auto mu = gauss->GetParameter(1);
auto sigma = gauss->GetParameter(2);
auto mu_fit_info = TString::Format("#mu = %.3f", mu);
auto su_fit_info = TString::Format("#sigma = %.3f", sigma);

legend->AddEntry(gauss, mu_fit_info.Data(), "l");
legend->AddEntry(pull_prt[name], su_fit_info.Data(), "");
}
}

if (predicted) {
std::string drawOptions = (smoothed or filtered) ? "same" : "";
auto drawOptions = (smoothed or filtered) ? "pe same" : "pe";

auto scale = 1. / pull_prt[name]->Integral("width");
pull_prt[name]->Scale(scale);
pull_prt[name]->Draw(drawOptions);
pull_prt[name]->Write();

pull_prt[vlID + paramNames.at(ipar)]->DrawNormalized(
drawOptions.c_str());
pull_prt[vlID + paramNames.at(ipar)]->Write();
legend->AddEntry(pull_prt[name], "predicted", "pe");

legend->AddEntry(pull_prt[vlID + paramNames.at(ipar)], "predicted",
"lp");
if (fitPredicted) {
pull_prt[name]->Fit("gaus", "q", "same");
TF1* gauss = pull_prt[name]->GetFunction("gaus");
gauss->SetLineColorAlpha(kRed, 0.5);

auto mu = gauss->GetParameter(1);
auto sigma = gauss->GetParameter(2);
auto mu_fit_info = TString::Format("#mu = %.3f", mu);
auto su_fit_info = TString::Format("#sigma = %.3f", sigma);

legend->AddEntry(gauss, mu_fit_info.Data(), "l");
legend->AddEntry(pull_prt[name], su_fit_info.Data(), "");
}
}

// Reference standard normal pdf
auto ref = new TF1("ref", "ROOT::Math::normal_pdf(x,1,0)", -5, 5);
ref->SetLineColor(kGreen);
ref->Draw("same");
legend->AddEntry(ref, "#mu = 0 #sigma = 1", "l");

legend->SetBorderSize(0);
legend->SetFillStyle(0);
legend->SetTextFont(42);
Expand Down