ROOT: Analyzing Petabytes of Data, Scientifically

借用 ROOT 官網上的介紹詞作標題。ROOT 是歐洲核子研究組織（CERN）開發的資料處理軟體，主要開發語言是 C++，使用時也是使用 C++ 以 include header 調用的。

我從大二就開始使用 Python 作為我做數據處理的工具， Python 有著方便編寫、不需編譯等優點，數據計算完後還可以直接用 Matplotlib 繪圖輸出向量圖檔，整個流程非常順暢。後來我在計算物理導論接觸到了 ROOT，當時我被 C++ 編譯器配置的繁瑣過程所困擾，並沒有成功利用這個強大的 Package。

最近因為一些偶然的機會再次接觸到用 C++ 做數據處理，花了一個下午終於完成了 ROOT 的環境配置；接下來我就介紹一下 ROOT 的優缺點吧。

速度快

ROOT 在使用時都以 C++ code 的形式調用的，在經過編譯後速度是 Python 的上萬倍，在處理大量資料、大量 fitting 時效果十分顯著。

將 C++ 轉換成免編譯語言

ROOT 的另一大優勢就是可以讓 C++ 成為不用編譯的語言，以腳本的方式執行，詳細請看以下例子：

Python 運作的方式：直譯器

ROOT 以直譯器的方式執行 C++ code

程式碼量較大

使用 C++ 的一大劣勢就是需要在 coding 的時間大幅度增加，加上 C++ 各種人類很難學明白的功能；同一個功能在 Python 跟 C++ 實現起來有著天壤之別。以下是用 Python + Numpy + Matplotlib + iminuit 跟使用 ROOT 實現同功能的程式碼比較。基本上是給定數據、做 curve fitting、繪圖。可以比較一下 coding 工作量跟成果：

import numpy as np
import matplotlib.pyplot as plt
from iminuit import Minuit

P = 84.35
B = 50
H = np.array([30.45, 24.8, 20.25, 19.35, 16])
V = np.array([20, 15, 10, 10, 10])
T = ([41.96, 45.13, 36.68, 41.54, 55.38])
Q = (V / 1000) / T
QT = (2 / 3) * np.sqrt(2 * 9.81) * (B / 1000) * (H / 1000)**(1.5)

vx = H
vy = Q
verr = np.ones(vx.size)

def func(K, n, x_ex):
    f = K * x_ex**n
    return f

def chi_sq(K, n):
    chi2 = (((func(K, n, vx) - vy) / verr) ** 2).sum()
    return chi2

m = Minuit(chi_sq, K = 3.0e-6, n = 1.5)
m.migrad()

K, n = m.values
chi2 = chi_sq(K, n)
span = vx[-1] - vx[0]
x_f = np.arange(vx[0] - 0.05 * span, vx[-1] + 0.05 * span, span / 8192)
y_f = func(K, n, x_f)
print(m.values, chi2)

font1 = {'family':'times','color':'black','size': 20}
font2 = {'family':'times','color':'black','size': 15}

Cd = K * (3 / 2) * (1 / 50) * (1 / np.sqrt(2 * 9.81)) * 10**(6 + n)

ABCDstring = "K = " + str(round(K, 10)) + ", n = " + str(round(n, 4)) + ", $C_d$ = " + str(round(Cd, 3))

fig, ax1 = plt.subplots(1)
fig.suptitle('')
fig.set_size_inches(8, 6)

ax1.plot(vx, vy, '+' , label = 'Experimental Value')

ax1.plot(x_f, y_f, '--', label = 'Regression Value \n' + ABCDstring)

ax1.set_xlabel("Water Surface Height (mm)", fontdict = font2)
ax1.set_ylabel("Discharge Quantity (m³/s)", fontdict = font2)


ax1.legend(loc = 0)
ax1.grid()
plt.savefig("Rectangular_2.pdf" , dpi = 900, bbox_inches = 'tight')
plt.show()

print(np.argmax(y_f))
print(x_f[np.argmax(y_f)])

print(Q)
print(QT)

#include "TGraphErrors.h"
#include "TCanvas.h"
#include "TF1.h"
#include "TAxis.h"
#include "TMinuit.h"
#include "Math/Minimizer.h"
#include "Math/Factory.h"
#include "Math/Functor.h"
#include "TLegend.h"
#include "TColor.h"
#include <iostream>

const int nPoints = 5;
double P = 84.35;
double B = 50;
double H[nPoints] = {30.45, 24.8, 20.25, 19.35, 16};
double V[nPoints] = {20, 15, 10, 10, 10};
double T[nPoints] = {41.96, 45.13, 36.68, 41.54, 55.38};
double Q[nPoints];
double QT[nPoints];
double verr[nPoints];


double func(const double *x, const double *p) {
    return p[0] * pow(x[0], p[1]);
}

double chi2(const double *par) {
    double sum = 0.0;
    for (int i = 0; i < nPoints; i++) {
        double delta = (func(&H[i], par) - Q[i]) / verr[i];
        sum += delta * delta;
    }
    return sum;
}

int main()
{
    for (int i = 0; i < nPoints; i++)
    {
        Q[i] = (V[i] / 1000) / T[i];
        QT[i] = (2.0 / 3.0) * sqrt(2.0 * 9.81) * (B / 1000.0) * pow((H[i] / 1000.0), 1.5);
        verr[i] = 1.0;
    }

    ROOT::Math::Minimizer *min = ROOT::Math::Factory::CreateMinimizer("Minuit2", "");
    
    ROOT::Math::Functor f(&chi2, 2);
    min->SetFunction(f);
    min->SetVariable(0, "K", 3.0e-6, 1.0e-7);
    min->SetVariable(1, "n", 1.5, 0.01);
    min->Minimize();

    const double *xs = min->X();
    std::cout << "K = " << xs[0] << ", n = " << xs[1] << std::endl;
    Int_t colorData = TColor::GetColor(57, 117, 179);
    Int_t colorFit = TColor::GetColor(242, 135, 33);

    TCanvas c1("c1", "", 800, 600);
    TGraphErrors gr(nPoints, H, Q, nullptr, nullptr);
    gr.SetMarkerStyle(2);
    gr.SetMarkerColor(colorData);
    gr.SetTitle(";Water Surface Height (mm);Discharge Quantity (m^3/s)");
    gr.Draw("AP");

    gr.GetXaxis()->SetTitle("Water Surface Height (mm)");
    gr.GetYaxis()->SetTitle("Discharge Quantity (m^{3}/s)");
    gr.GetXaxis()->CenterTitle();
    gr.GetYaxis()->CenterTitle();
    gr.GetXaxis()->SetTitleFont(132);
    gr.GetYaxis()->SetTitleFont(132);
    
    gPad->SetGridx(1);  // Enable grid for X axis
    gPad->SetGridy(1);  // Enable grid for Y axis
    gPad->Update();

    TF1 fitFunction("fitFunction", func, H[0], H[nPoints-1], 2);
    fitFunction.SetParameters(xs[0], xs[1]);
    fitFunction.SetLineColor(colorFit);
    fitFunction.SetLineStyle(2);
    fitFunction.Draw("SAME");


    // Leagend
    double Cd = xs[0] * (3.0 / 2.0) * (1.0 / 50.0) * (1.0 / sqrt(2.0 * 9.81)) * pow(10.0, 6.0 + xs[1]);
    TString ABCDstring = TString::Format("K = %.4e, n = %.4f, C_{d} = %.3f", xs[0], xs[1], Cd);
    
    TLegend *legend = new TLegend(0.4, 0.125, 0.875, 0.25); // x1, y1, x2, y2
    legend->SetFillColorAlpha(0, 0.8); // kGray for gray color, 0.5 for 50% transparency
    legend->AddEntry(&gr, "Experimental Value", "P");
    
    legend->AddEntry(&fitFunction, "Regression Value", "L");
    legend->AddEntry((TObject*)0, ABCDstring.Data(), ""); // Display the regression formula
    
    legend->Draw();

    c1.SaveAs("./Rectangular_2_2.pdf");

    return 0;
    
}

用 Matplotlib 畫出來的結果

用 ROOT 畫出來的結果

可以看出算出來的結果是差不多的，同時是實際執行時 ROOT 的速度快了好幾倍，但可以看出繪製出來的圖片略顯粗糙；可以說是各有優缺點。