Advanced Computing Platform for Theoretical Physics

Skip to content

GitLab

Commit 838a8da6 authored by Yu-Chen Ding's avatar Yu-Chen Ding
Browse files

ratio plot

parent 41b0fa56
Hide whitespace changes
Inline Side-by-side
include/crdb.h
View file @ 838a8da6
......@@ -148,7 +148,11 @@ class CRPlot {
bool log_x, log_y;
bool grid_x, grid_y;
std::string ratio_ref;
double ratio_min, ratio_max;
int canvas_width, canvas_height;
double divide_ratio;
// todo: normalized comparisons
......
src/crplot.cc
View file @ 838a8da6
#include "crutil.hh"
#include <crdb.h>
#include <iostream>
#include <array>
#include <cmath>
......@@ -9,9 +10,13 @@
#include <TColor.h>
#include <TFile.h>
#include <TGaxis.h>
#include <TGraphAsymmErrors.h>
#include <TH1D.h>
#include <TLegend.h>
#include <TPad.h>
#include "loginterp"
namespace _CRPLOT {
template <size_t N>
......@@ -75,8 +80,10 @@ CRPlot::CRPlot(const std::string& title, int canv_w, int canv_h)
log_x(true), log_y(true), //
grid_x(true), grid_y(true), //
scale_index(0), legend_phase(4), //
ratio_min(NAN), ratio_max(NAN), //
canvas_width(canv_w), //
canvas_height(canv_h) {}
canvas_height(canv_h), //
divide_ratio(0.3) {}
using Style = CRPlot::Style;
void CRPlot::AddData(const CRData& data, const Style& style,
......@@ -196,48 +203,140 @@ double AdjustYMax(const CRPlot& plot) {
}
std::shared_ptr<TCanvas> CRPlot::Plot(void) const {
double x_min_plot = std::isnan(x_min) ? AdjustXMin(*this) : x_min;
double x_max_plot = std::isnan(x_max) ? AdjustXMax(*this) : x_max;
double y_min_plot = std::isnan(y_min) ? AdjustYMin(*this) : y_min;
double y_max_plot = std::isnan(y_max) ? AdjustYMax(*this) : y_max;
const bool with_ratio = !ratio_ref.empty();
double x_min_plot = std::isnan(x_min) ? AdjustXMin(*this) : x_min;
double x_max_plot = std::isnan(x_max) ? AdjustXMax(*this) : x_max;
double y_min_plot = std::isnan(y_min) ? AdjustYMin(*this) : y_min;
double y_max_plot = std::isnan(y_max) ? AdjustYMax(*this) : y_max;
double y_min_ratio = with_ratio && std::isnan(ratio_min) ? 0.5 : ratio_min;
double y_max_ratio =
with_ratio && std::isnan(ratio_max) ? (log_y ? 2.0 : 1.5) : ratio_max;
if (std::isnan(x_min_plot) || std::isnan(x_max_plot) ||
std::isnan(y_min_plot) || std::isnan(y_max_plot)) {
throw std::runtime_error("invalid plot range");
}
auto pcanv = std::make_shared<TCanvas>("canv", "canv", canvas_width, canvas_height);
// calclulate canvas layout
const double label_size = canvas_width * 0.02;
const double title_size = canvas_width * 0.03;
const int label_font = 133;
const int title_font = 133;
// overall canvas
auto pcanv =
std::make_shared<TCanvas>("canv", "canv", canvas_width, canvas_height);
auto& canv = *pcanv;
canv.SetFillColor(0);
// canv.SetFillColor(0);
// canv.SetGrid(grid_x, grid_y);
canv.SetBorderMode(0);
canv.SetLogx(log_x);
canv.SetLogy(log_y);
canv.SetGrid(grid_x, grid_y);
canv.SetTickx(1);
canv.SetTicky(1);
canv.SetFrameBorderMode(0);
// TH1D hframe("hframe", title.c_str(), 1, x_min_plot, x_max_plot);
TH1D* phframe =
new TH1D("hframe", title.c_str(), 1, x_min_plot, x_max_plot);
phframe->SetBit(kCanDelete); // prevent leak
auto& hframe = *phframe;
hframe.SetMinimum(y_min_plot);
hframe.SetMaximum(y_max_plot);
hframe.SetStats(0);
hframe.SetXTitle(x_label.c_str());
hframe.SetYTitle(y_label.c_str());
hframe.Draw();
// th1d for set plot range range
TH1D *phframe = nullptr, *phframe_ratio = nullptr;
phframe = new TH1D("hframe", title.c_str(), 1, x_min_plot, x_max_plot);
phframe->SetMinimum(y_min_plot);
phframe->SetMaximum(y_max_plot);
if (with_ratio) {
phframe_ratio = new TH1D("hframe_ratio", "", 1, x_min_plot, x_max_plot);
phframe_ratio->SetMinimum(y_min_ratio);
phframe_ratio->SetMaximum(y_max_ratio);
}
for (auto phf : {phframe, phframe_ratio}) {
if (!phf) {
continue;
}
phf->SetBit(kCanDelete); // prevent leak
phf->SetStats(0);
phf->GetYaxis()->SetDecimals(!log_y);
phf->GetYaxis()->SetLabelFont(label_font);
phf->GetYaxis()->SetLabelSize(label_size);
}
phframe->GetYaxis()->SetTitleFont(title_font);
phframe->GetYaxis()->SetTitleSize(title_size);
phframe->GetYaxis()->SetTitle(y_label.c_str());
auto xax = (phframe_ratio ? phframe_ratio : phframe)->GetXaxis();
xax->SetTitleSize(title_size);
xax->SetTitleFont(title_font);
xax->SetTitle(x_label.c_str());
xax->SetLabelFont(label_font);
xax->SetLabelSize(label_size);
xax->SetTitleOffset(3.);
xax->SetDecimals(!log_x);
TPad *pad_plot = nullptr, *pad_ratio = nullptr;
yuc::log_interpolator interp_fcn;
if (with_ratio) {
pad_plot = new TPad("pad_plot", "pad_plot", 0, divide_ratio, 1, 1);
pad_ratio = new TPad("pad_ratio", "pad_ratio", 0, 0, 1, divide_ratio);
pad_plot->SetBottomMargin(0);
pad_ratio->SetTopMargin(0);
pad_ratio->SetBottomMargin(0.1 / divide_ratio);
// prepare interpolation function TODO: linear interp for linear data
auto ref_graph = _graphs.find(ratio_ref);
auto ref_egraph = _egraphs.find(ratio_ref);
auto ref_band = _bands.find(ratio_ref);
if (ref_graph != _graphs.end()) {
interp_fcn.insert(ref_graph->second.first.x,
ref_graph->second.first.y);
} else if (ref_egraph != _egraphs.end()) {
interp_fcn.insert(ref_egraph->second.first.x(),
ref_egraph->second.first.y());
} else if (ref_band != _bands.end()) {
interp_fcn.insert(ref_band->second.first.x,
ref_band->second.first.y);
} else {
throw std::runtime_error("ratio reference not found: '" +
ratio_ref + "'");
}
} else {
pad_plot = new TPad("pad_plot", "pad_plot", 0, 0, 1, 1);
}
for (auto spad : {pad_plot, pad_ratio}) {
if (!spad) {
continue;
}
spad->SetLogx(log_x);
spad->SetLogy(log_y);
spad->SetGrid(grid_x, grid_y);
spad->SetFillColor(0);
pcanv->cd();
spad->Draw();
}
pad_plot->cd();
phframe->Draw();
if (with_ratio) {
pad_ratio->cd();
phframe_ratio->Draw();
}
// prepare legends
size_t n_entry = 0;
auto leg = new TLegend(0.1, with_ratio ? 0 : 0.1, 0.9, 0.9);
leg->SetBit(kCanDelete); // prevent leak
// do not delete TGraph and TGraphAsymmErrors until finished!
std::vector<TGraph*> graphs;
std::vector<TGraphAsymmErrors*> egraphs;
std::vector<TGraph*> ratio_graphs;
std::vector<TGraphAsymmErrors*> ratio_egraphs;
for (const auto& d : _bands) {
const auto& band = d.second.first;
const auto style = d.second.second;
const size_t n = band.x.size();
auto g = new TGraphAsymmErrors(n);
g->SetBit(kCanDelete); // prevent leak
const auto& band = d.second.first;
const auto style = d.second.second;
const std::string draw_flag = style.size > 0 ? "EL3" : "E3";
const std::string draw_leg_flag = style.size > 0 ? "fl" : "f";
const size_t n = band.x.size();
auto g = new TGraphAsymmErrors(n);
auto rg = with_ratio ? new TGraphAsymmErrors(n) : nullptr;
g->SetTitle(d.first.c_str());
for (size_t i = 0; i < n; ++i) {
double x = band.x[i];
......@@ -256,55 +355,98 @@ std::shared_ptr<TCanvas> CRPlot::Plot(void) const {
g->SetPoint(i, x, y * fac);
g->SetPointError(i, 0, 0, el * fac, eu * fac);
if (rg) {
double ref_y = interp_fcn(x);
rg->SetPoint(i, x, y / ref_y);
rg->SetPointError(i, 0, 0, el / ref_y, eu / ref_y);
}
}
if (style.fill_style / 1000 == 1) {
double alpha = (style.fill_style % 1000 % 101) * .01;
g->SetFillColorAlpha(
style.fill_color < 0 ? style.color : style.fill_color, alpha);
g->SetFillStyle(1001);
} else {
g->SetFillColor(style.fill_color < 0 ? style.color
: style.fill_color);
g->SetFillStyle(style.fill_style);
}
if (style.size > 0) {
// draw center line
g->SetLineWidth(style.size);
g->SetLineColor(style.color);
g->SetLineStyle(style.style);
g->Draw("same CE3");
} else {
g->SetLineWidth(0);
g->Draw("same E3");
for (auto sg : {g, rg}) {
if (!sg) {
continue;
}
sg->SetBit(kCanDelete); // prevent leak
if (style.fill_style / 1000 == 1) {
double alpha = (style.fill_style % 1000 % 101) * .01;
sg->SetFillColorAlpha(style.fill_color < 0 ? style.color
: style.fill_color,
alpha);
sg->SetFillStyle(1001);
} else {
sg->SetFillColor(style.fill_color < 0 ? style.color
: style.fill_color);
sg->SetFillStyle(style.fill_style);
}
if (style.size > 0) {
// draw center line
sg->SetLineWidth(style.size);
sg->SetLineColor(style.color);
sg->SetLineStyle(style.style);
} else {
sg->SetLineWidth(0);
sg->SetLineColor(0);
sg->SetLineStyle(0);
}
}
pad_plot->cd();
g->Draw(("same " + draw_flag).c_str());
egraphs.push_back(g);
if (rg) {
pad_ratio->cd();
rg->Draw(("same " + draw_flag).c_str());
ratio_egraphs.push_back(g);
}
leg->AddEntry(g, d.first.c_str(), draw_leg_flag.c_str());
++n_entry;
}
// plot curves
for (const auto& d : _graphs) {
const auto& curve = d.second.first;
const auto style = d.second.second;
auto g = new TGraph(curve.x.size());
g->SetBit(kCanDelete); // prevent leak
const std::string draw_flag = "l";
const auto& curve = d.second.first;
const auto style = d.second.second;
const size_t n = curve.x.size();
auto g = new TGraph(curve.x.size());
auto rg = with_ratio ? new TGraph(n) : nullptr;
g->SetTitle(d.first.c_str());
for (size_t i = 0; i < curve.x.size(); ++i) {
double x = curve.x[i], y = curve.y[i],
fac = std::pow(x, scale_index);
g->SetPoint(i, x, y * fac);
if (rg) {
double ref_y = interp_fcn(x);
rg->SetPoint(i, x, y / ref_y);
}
}
g->SetLineColor(style.color);
g->SetLineStyle(style.style);
g->SetLineWidth(style.size < 0 ? 2 : style.size);
g->Draw("same l");
for (auto sg : {g, rg}) {
if (!sg) {
continue;
}
sg->SetBit(kCanDelete); // prevent leak
sg->SetLineColor(style.color);
sg->SetLineStyle(style.style);
sg->SetLineWidth(style.size < 0 ? 2 : style.size);
}
pad_plot->cd();
g->Draw(("same " + draw_flag).c_str());
graphs.push_back(g);
if (rg) {
pad_ratio->cd();
rg->Draw(("same " + draw_flag).c_str());
ratio_graphs.push_back(g);
}
leg->AddEntry(g, d.first.c_str(), draw_flag.c_str());
++n_entry;
}
// plot data
for (const auto& d : _egraphs) {
const auto data = d.second.first;
const auto style = d.second.second;
auto g = new TGraphAsymmErrors(data.size());
g->SetBit(kCanDelete); // prevent leak
const std::string draw_flag = "ep";
const auto data = d.second.first;
const auto style = d.second.second;
const size_t n = data.size();
auto g = new TGraphAsymmErrors(data.size());
auto rg = with_ratio ? new TGraphAsymmErrors(n) : nullptr;
g->SetTitle(d.first.c_str());
for (size_t i = 0; i < data.size(); ++i) {
auto dp = data[i];
......@@ -317,35 +459,63 @@ std::shared_ptr<TCanvas> CRPlot::Plot(void) const {
dp.x - dp.x_lower, //
dp.x_upper - dp.x, //
dp.ey1 * fac, dp.ey2 * fac);
if (rg) {
double ref_y = interp_fcn(dp.x);
rg->SetPoint(i, dp.x, dp.y / ref_y);
rg->SetPointError(i, dp.x - dp.x_lower, dp.x_upper - dp.x,
dp.ey1 / ref_y, dp.ey2 / ref_y);
}
}
for (auto sg : {g, rg}) {
if (!sg) {
continue;
}
sg->SetBit(kCanDelete); // prevent leak
sg->SetFillColor(0);
sg->SetLineColor(style.color);
sg->SetMarkerColor(style.color);
sg->SetMarkerStyle(style.style);
sg->SetLineWidth(style.size < 2 ? 1 : (style.size - 1));
sg->SetMarkerSize(style.size < 0 ? 2 : style.size);
}
g->SetFillColor(0);
g->SetLineColor(style.color);
g->SetMarkerColor(style.color);
g->SetMarkerStyle(style.style);
g->SetLineWidth(style.size < 2 ? 1 : (style.size - 1));
g->SetMarkerSize(style.size < 0 ? 2 : style.size);
g->Draw("same ep");
pad_plot->cd();
g->Draw(("same " + draw_flag).c_str());
egraphs.push_back(g);
if (rg) {
pad_ratio->cd();
rg->Draw(("same " + draw_flag).c_str());
ratio_egraphs.push_back(g);
}
leg->AddEntry(g, d.first.c_str(), draw_flag.c_str());
++n_entry;
}
size_t nLedgend = egraphs.size() + graphs.size();
TLegend* l;
switch (legend_phase) {
case 1:
l = canv.BuildLegend(0.6, 0.9 - nLedgend * .03, 0.9, 0.9);
// l = pad_plot->BuildLegend(0.6, 0.9 - nLedgend * .03, 0.9, 0.9);
leg->SetX1(0.6);
leg->SetY1(leg->GetY2() - n_entry * 0.03);
break;
case 2:
l = canv.BuildLegend(0.1, 0.9 - nLedgend * .03, 0.4, 0.9);
// l = pad_plot->BuildLegend(0.1, 0.9 - nLedgend * .03, 0.4, 0.9);
leg->SetX2(0.4);
leg->SetY1(leg->GetY2() - n_entry * 0.03);
break;
case 3:
l = canv.BuildLegend(0.1, 0.1, 0.4, 0.1 + nLedgend * .03);
// l = pad_plot->BuildLegend(0.1, 0.1, 0.4, 0.1 + nLedgend * .03);
leg->SetX2(0.4);
leg->SetY2(leg->GetY1() + n_entry * 0.03);
break;
case 4:
l = canv.BuildLegend(0.6, 0.1, 0.9, 0.1 + nLedgend * .03);
// l = pad_plot->BuildLegend(0.6, 0.1, 0.9, 0.1 + nLedgend * .03);
leg->SetX1(0.6);
leg->SetY2(leg->GetY1() + n_entry * 0.03);
break;
}
l->GetListOfPrimitives()->RemoveAt(0);
l->SetBit(kCanDelete); // prevent leak
// TODO multi-column legend
pad_plot->cd();
leg->Draw();
return pcanv;
......
src/loginterp 0 → 100644
View file @ 838a8da6
// version 2.0
#pragma once
#include <cmath>
#include <set>
namespace yuc {
class log_interpolator {
protected:
struct data_point {
double x, y;
bool operator<(const data_point& d) const { return x < d.x; }
};
protected:
std::set<data_point> data;
public:
log_interpolator(void) {}
template <typename V> log_interpolator(const V& xa, const V& ya) {
insert(xa, ya);
}
public:
void clear(void) { data.clear(); }
void insert(double x, double y) {
if (x >= 0 && y >= 0) {
data.insert({std::log(x), std::log(y)});
}
}
template <typename V> void insert(const V& xa, const V& ya) {
for (size_t i = 0; i < xa.size(); ++i) {
this->insert(xa[i], ya[i]);
}
}
public:
double operator()(double x) const {
if (x < 0 || data.size() < 2) {
return NAN;
}
x = std::log(x);
auto it1 = data.upper_bound({x}), it2 = std::prev(it1);
if (it1 == data.begin()) {
it2 = std::next(it1);
} else if (it1 == data.end()) {
it1 = std::prev(it2);
}
const double x1 = it1->x, y1 = it1->y;
const double x2 = it2->x, y2 = it2->y;
if (x1 == x2) {
return std::exp(0.5 * (y1 + y2));
}
if (y1 == -HUGE_VAL || y2 == -HUGE_VAL) {
if (y1 == y2) {
return 0.;
} else {
// use log-linear interpolation if one y is zero
return std::exp(y1) +
(x - x1) * (std::exp(y2) - std::exp(y1)) / (x2 - x1);
}
}
return std::exp(y1 + (x - x1) * (y2 - y1) / (x2 - x1));
// (power-law extrapolating behavior)
}
template <typename V> V operator()(const V& xs) const {
V ys(xs.size());
for (size_t i = 0; i < xs.size(); ++i) {
ys[i] = this->operator()(xs[i]);
}
return ys;
}
};
};
// vi:ft=cpp
test/CMakeLists.txt
View file @ 838a8da6
......@@ -17,6 +17,12 @@ add_test(
NAME test:plot
COMMAND test_plot "${Test_DATA}/BC_EKN.usine" "${Test_DATA}/fitted_BC.tsv" "${Test_DATA}/fitted_BC_band.tsv")
add_executable(test_plot_with_ratio "test_plot_with_ratio.cc")
target_link_libraries(test_plot_with_ratio crdb ${ROOT_LIBRARIES})
add_test(
NAME test:plot_with_ratio
COMMAND test_plot_with_ratio "${Test_DATA}/BC_EKN.usine" "${Test_DATA}/fitted_BC.tsv" "${Test_DATA}/fitted_BC_band.tsv")
add_executable(test_fetch "test_fetch.cc")
target_link_libraries(test_fetch crdb ${ROOT_LIBRARIES})
add_test(NAME test:fetch_B_C
......
test/test_plot_with_ratio.cc 0 → 100644
View file @ 838a8da6
#include <iostream>
#include <crdb.h>
int main(int argc, char* argv[]) {
if (argc != 4) {
throw std::runtime_error(std::string(argv[0]) + " <data file> <curve file> <band file>");
}
CRPlot plot("B/C Test Plot");
CRDataset ds;
ds.parse_file(argv[1]);
auto data = ds.find_exact("B/C", "AMS02");
plot.AddData(*data, {2, 20, 1}, "AMS02");
// plot.AddCurve(argv[2], {3, 2, 2}, "fitted B/C");
plot.AddBand(argv[3], {4, 2, 2, 1020}, "fitted B/C");
plot.log_y = false;
plot.x_label = "E_{kin}/nuc [GeV/n]";
plot.y_label = "B/C";
plot.grid_x = plot.grid_y = true;
plot.legend_phase = 1;
plot.ratio_ref = "fitted B/C";
auto canvas = plot.Plot();
canvas->SaveAs("test_plot_with_ratio.pdf");
canvas->SaveAs("test_plot_with_ratio.png");
canvas->SaveAs("test_plot_with_ratio.C");
return 0;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Provided by Institute of Theoretical Physics, Chinese Academy of Sciences. Please contact support@itp.ac.cn for more information.