用C++做資料分析 | CERN ROOT 教學[03]

簡介

在上一篇文章中我們介紹了基本資料分析的操作方法。在這篇文章中我們會詳細介紹在資料分析中常見的一個功能：函數擬合。

函數

ROOT的TF1可以用來表示任意的一元函數。他的自變數的名稱永遠為x。而且你可以在裡面使用任意的C/C++函數(包括ROOT自帶的TMath函數庫)。

例如創造一個範圍在[-10, 10]之間的sinc函數。

root [0] f1 = new TF1("f1", "sin(x)/x", -10, 10);
root [1] f1->Draw();


sin(x)/x

因為你可以用任何 C/C++的函數。所以也可以自己定義一個std::function或是lambda function後使用她。值得注意的是因為C++17支援Class template type deduction而且ROOT自動匯入std的命名空間。所以原本在C++11宣告一個std::function的複雜語句auto my_func = std::function<double(double)>([](double){...})可以簡化成my_func = function([](double){...})。

root [2] my_sinc = function([](double x){ return sin(x)/x; });
root [3] f1 = new TF1("f1", "my_sinc(x)", -10, 10);
root [4] f1->Draw();


sin(x)/x

說到這裡熟悉C++的讀者應該會想到：那具有function call operator的類別是否也可以當作函數使用呢？答案是可以的。我在這裡就不多做說明，不過應該很容易成功。

有了函數，還需要能夠設定函數的各種參數。像是想要表達一個任意的線性函數y = ax+b的時候需要能夠指定固定的a與b。我們可以在函數的表達式中加入[0]來代表第0個參數、[1]代表第一個，以此類推。並在繪製圖形前用TF1::SetParameters來設定數值(其型別固定為double)

root [5] f1 = new TF1("f1", "[0]*x+[1]", -10, 10);
root [6] f1->SetParameters(-1, -5);
root [7] f1->Draw();


y = ax+b，a = -1, b = -5

你也可以使用變數名稱來指定變數。這種方法比較好寫也比較容易除錯。但不像使用索引一樣使用SetParameters一次指定一堆變數。

root [8] f1 = new TF1("f1", "exp(x+[a])", -10, 10);
root [9] f1->SetParameter("a", 2);
root [10] f1->Draw();


exp(x+a)

當需要取得變數的數值的時候可以使用GetParameter來取得函數。

root [11] f1->GetParameter("a")
(double) 2.0000000

函數擬合

講完函數之後終於可以討論如何使用函數擬合了！首先。我們先生出一個常態分佈。

root [0] h1 = new TH1F("h1", "histogram", 60, -3, 3);
root [1] h1->FillRandom("gaus", 100000);
root [2] h1->Draw();


常態分佈

並且如前一篇文章所述，使用h1->Fit("gaus");將高斯分佈函數擬合上去。

root [4] h1->Fit("gaus");
 FCN=58.2284 FROM MIGRAD    STATUS=CONVERGED      62 CALLS          63 TOTAL
                     EDM=8.58021e-13    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  Constant     4.00036e+03   1.57642e+01   4.75510e-02  -9.85101e-09
   2  Mean        -1.17823e-03   3.20398e-03   1.20393e-05  -2.62278e-05
   3  Sigma        9.99354e-01   2.41238e-03   2.44449e-06  -1.71898e-03


將高斯分佈擬合至高斯分佈上

之後我們可以使用GetFunction取得擬合後函數本身並取得函數的所有參數。在此函數的名稱為你所使用的分佈，而參數的名稱在上面擬合後的輸出結果中。

root [5] f1 = h1->GetFunction("gaus");
root [6] f1->GetParameter("Constant")
(double) 4000.3641
root [7] f1->GetParameter("Sigma")
(double) 0.99935417

擬合自定義函數

當然我們不只能擬合到預先定義好的函數上，也可以擬合到自定義的函數。不過有兩個前提：要擬合的函數必須帶有參數，且函數必須是determinstic的。所以ROOT可以幫你擬合[scale]*x+[bias]+exp(x*[exp_scale])但沒辦法擬合(rand()%1000)*x*[scale]。

要擬合自定義函數需要先定義一個TF1函數。然後呼叫TH1::Fit將函數擬合上去。

root [7] f2 = new TF1("f2", "[scale]*sin(x)/x+[bias]");
root [8] h1->Fit(f2);
 FCN=14310 FROM MIGRAD    STATUS=CONVERGED      30 CALLS          31 TOTAL
                     EDM=7.17638e-19    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  bias        -3.32813e+02   3.54589e+00   1.38923e-01   2.29136e-10
   2  scale        2.85750e+03   1.11800e+01   4.38015e-01  -4.15280e-11