Root - framework (nejen) pro analýzu dat

23.12.2007 20:59 | Přečteno: 1068×

Krátký úvod do jednoho objektově orientovaného frameworku pro analýzu dat.

Právě dělám jeden prográmek, ve kterém využívám Root. Chtěl jsem se podělit o některé zkušenosti. Jelikož ale Root nepatří k nejznámějším, prvně se zde pokusím popsat oč se vůbec jedná a ukázat jednoduchý příklad použití.

Jak již bylo uvedeno výše, Root je objektově orientovaný, v C++ napsaný, framework pro analýzu dat. Je primárně určený pro analýzu dat ve fyzice vysokých energií, což ale nijak nebrání tomu použít některé jeho funkce i pro jiná data. Vyvíjen je v CERNu, kde by také měl být používán na zpracování dat z LHC (což mimo jiné znamená, že vývoj Rootu a jeho podpora v dohledné době určitě neskončí). Vydáván je pod GNU/GPL licencí. Existují interface Rootu pro Python - pyroot a údajně i pro Ruby.

Root se skládá ze dvou hlavních částí. Jednak ze samotného Rootu a druhak z CINTu - intepretru C++ rozšířeného o Root. Příkazem root spustíme interpretr. Ukončíme jej příkazem .q. Makro v externím souboru spustíme příkazem .x nazev_souboru. Pro větší programy makra nejsou příliš vhodná, ale na úvodní osahaní Rootu zcela dostatečná. Pro to ale stačí i samotná příkazová řádka interpretru.

Nejjednodušší věc, s kterou můžeme začít, jsou jednorozměrná data, která zobrazíme jako histogram. Prvně vytvoříme jednorozměrný histogram (příkazy buď můžeme rovnou psát do příkazové řádky interpretru, do souboru, který pak spustíte pomocí .x soubor nebo jako program, který pak zkompilujete, viz závěr):
TH1F *h = new TH1F("h","prvni histogram",100,100.,250.); Třída TH1F slouží k vytvoření jednorozměrného histogramu, který bude obsahovat data typu float. První parametr konstruktoru je název histogramu, druhý jeho popis, který pak můžeme vykreslit společně s histogramem. Třetí patametr udává počet binů histogramu, předposlední parametr dolní mez a poslední mez horní (šířka binu pak tedy je (horniMez - dolniMez)/pocetBinu). Máme tedy jednorozměrný histogram, který prvně musíme naplnit daty. Prvně je ale někde musíme vzít. Pokud zrovna žádná nemáme, tak si je můžeme třeba vygenerovat (Root nabízí hned několik generátorů (pseudo)náhodných čísel):
Trandom *r = new TRandom();
for(int i=0;i<10000;i++){
  h->Fill(r->Gaus(170,15));
}
Nyní jsme tedy vygenerovali 10000 náhodných čísel podle gaussovského rozdělení se střední hodnotou 170 a rozptylem 15 (r->Gaus(170,15)). Těmito čísly jsem pak naplnili histogram (h->Fill(cislo);). Histogram můžeme zobrazit zavoláním metody Draw(): h->Draw(); V pravém rohu se nám při standardním nastavení zobrazí počet záznamů, průměr a RMS (v našem případě překvapivě něco jako 169,8 a 15,1). Co se vhledu týče, můžeme nastavit téměř cokoli a to buď přímo z našeho programu nebo naklikat díky editoru, který otevřeme z menu (View->Editor).

Vygenerujme nyní o trošku složitější data a pokusme se je nafitovat:
for(int i=0;i<10000;i++){
  h->Fill(r->Gaus(175,15));
  h->Fill(r->Gaus(140,10));
}
h->Draw();
I kdyby jsme o datech nic nevěděli, první, co by nás asi napadlo, by bylo zkusit tyto data nafitovat složením dvou gausovek. Zadefinujeme si tedy tuto funkci, hranatá závorka označuje n-tý parametr funkce, poslední dvě čísla zadávají interval, na kterém je funkce definována:
TF1* fc1 = new TF1("fc1","(1/(sqrt(2*TMath::Pi())*[0]))*[2]*exp(-(x-[1])**2/(2*[0]**2)) + (1/(sqrt(2*TMath::Pi())*[3]))*\ [5]* exp(-(x-[4])**2/(2*[3]**2))",100,250); Aby fit dopadnul rozumě, musíme aspoň některým parametrům nastavit počáteční hodnoty. Přibližnou hodnotu parametrů 1 a 3 snadno vidíme, 4 a 5 nastavíme na nějakou nenulovou hodnotu, např takto:
fc1->SetParameter(1, 140);
fc1->SetParameter(2,100);
fc1->SetParameter(4, 180);
fc1->SetParameter(5,100);
Fit provedeme jednoduše:
h->Fit("fc1"); Root provede fit a vypíše parametry, jejich chybu etc. a zobrazí fit. Parametry můžeme získat zavoláním metody GetParameters(double* params), s funkcí můžeme rovněž dál pracovat, např. můžeme vykreslit její integrál.

Celý program, který následně zkompilujeme, by mohl vypadat takto (plátno ještě rozdělíme a na konci uložíme jako gif):
#include <TApplication.h>
#include <TH1F.h>
#include <TF1.h>
#include <TCanvas.h>
#include <TRandom3.h>
#include <TMath.h>
int main(int argc, char **argv) {
   TApplication theApp("App", &argc, argv);
   TCanvas *c1 = new TCanvas("c1", "c1",15,45,699,499);
   c1->Divide(2,2);
  
   TH1F *h = new TH1F("h","prvni histogram",100,100.,250.);
   TRandom *r = new TRandom();
   for(int i=0;i<10000;i++){
   h->Fill(r->Gaus(175,15));
   h->Fill(r->Gaus(140,10));
   }
   c1->cd(1);
   h->Draw();
  
   TF1* fc1 = new TF1("fc1","(1/(sqrt(2*TMath::Pi())*[0]))*[2]*exp(-(x-[1])**2/(2*[0]**2)) + (1/(sqrt(2*TMath::Pi())*[3]))*[5]* exp(-(x-[4])**2/(2*[3]**2))",100,250);
  
   fc1->SetParameter(1, 140);
   fc1->SetParameter(2,100);
   fc1->SetParameter(4, 180);
   fc1->SetParameter(5,100);
  
   c1->cd(2);
   TH1F *hTmp = (TH1F*)h->Clone();
   hTmp->Draw();
   hTmp->Fit("fc1");
   c1->cd(3);
   fc1->Draw();
   c1->cd(4);
   fc1->DrawIntegral();
   c1->SaveAs("example1.gif");
  
   theApp.Run();
   return 0;
}
Program můžeme zkompilovat takto:
g++ -I `root-config --incdir` -o test example1.C `root-config --libs`
Program spustíme normálně, ale ukončit jej musíme buď z shellu pomocí Ctrl-c nebo z menu okna File->Quit Root. Jelikož jsme neimplementovali žádné naslouchání událostem, tak klasickým zavřením okna se program neukončí.

       

Anketa

Obrázky

Komentáře

[frdrx]

Ha, tohle by mělo jít dobře použít na analýzu píčků z cytometru, co? Musím na to kouknout.