CERN RooStat/RooFit

Python Roo系

CERN ROOT ライブラリを使用して、統計処理を「簡単に」行うためのはずのパッケージ群。実際には使いにくい（慣れてもあまり使いたくない）です。以下のチュートリアルでは

import ROOT as R

が宣言されていることを念頭においています。

基本的なクラス一覧

素人が作ったんかと思うくらいウンザリする。

クラス	説明	感想
RooPlot	何かを描画するためのフレーム。よくframeという変数に格納される	RooFrameとかの名前にしろよ
RooHist	TGraphAsymErrorsを継承していて、ビニングされたプロットを表現する	TGraphは基本ヒストグラムではないやろ
RooCurve	TGraphを継承していて、1次元グラフを表現する	カーブて...色々あるやろ...

RooFit

RooHist

hist.getFitRangeNEvt() # イベント数をprintする

基本思想

数式	RooFitのクラス	感想
$x$	RooRealVar
$\vec{x}$	RooArgSet
$f(x)$	RooAbsReal
$F(x)$	RooAbsPdf
$\int_{a}^{b} f(x)~dx$	RooRealIntegral

乱数をふる

自分で確率密度関数を定義して、それに従った乱数を振ることができる。ちなみに基本的に何でもチュートリアルは$ROOTSYS/tutorialsにあるので随時参照。ここでは$ROOTSYS/tutorials/roofit/のコードを改造したものを備忘録としてまとめる。

RooWorkspace::factory()関数を使う（おすすめ）

using namespace RooFit;
using namespace RooStats;

void generate_gaussian()
{
    RooWorkspace ws("ws");
    ws.factory("Gaussian::gaus(x[-10,10], mean[0], sigma[1])");

    RooDataSet* data = ws.pdf("gaus")->generate(*ws.var("x"), 1000);
    RooPlot *frame = ws.var("x")->frame();
    data->plotOn(frame);
    ws.pdf("model")->plotOn(frame);
    frame->Draw();
}

[解説]
まず、何よりも先にRooWorkspaceを作る。そして、乱数が従う確率密度関数を定義する（ここではガウス分布）。factoryでGaussian::<ユーザー定義の名前>という文法でガウス分布を作っていて、

G(\mu=0,\sigma=1) = \frac{1}{\sqrt{2\pi}\sigma}~e^{\left( \frac{x-\mu}{\sigma} \right)^2}

の関数を定義している。Gaussian::gausの引数は

• $x$：観測量。区間を定義することができ、ここでは$-10<x<10$の範囲で観測される場合を表現している
• mean：$\mu$。mean[1, -1, 1]みたいな書き方だったら、meanを1に固定して、その定数の範囲（プロットするときに使用する下限上限だと思う）を定義している意味になる。mean[-1,1]としたら$-1<\mu<1$の範囲の変数を意味することとなる。

RooWorkspaceを使わない

ちなみに全く同じことが以下のコードでもできる。やり方が何通りかあって、調べる資料によってバラバラなところがうっとおしい。

void generate_gaussian_2()
{
    RooRealVar x("x","x",-10,10);
    RooRealVar mean("mean","mean",0);
    RooRealVar sigma("sigma","sigma",1);
    RooGaussian gaus("gaus","gaus",x,mean,sigma);

    RooDataSet* data = gaus.generate(x, 1000);
    RooPlot* frame = x.frame();
    data->plotOn(frame);
    frame->Draw();
}

これを見ると、RooWorkspaceを使ったほうが良さそうな気がするのは確か。直接自分で複数のPDFを記述しようとするとRooRealVarが大量に増えていって、それらを管理するための自作クラスを導入したい気持ちに駆られると思うので。

パラメーター推定

興味のあるパラメーター（parameter of interest; POI）の尤もらしい値を推定する際に、Likelihoodを使用する。Likelihoodが最大となるような時に、そのPOIは尤もらしい値となる。また Maximum Likelihood Estimation (MLE) ではなく、Negative Log-likelihood が最小となる場合と解釈し直して計算することが多い（計算しやすい）。

Likelihood

確率密度関数（probability distribution function; pdf）をかけ合わせた確率”のようなもの”。

L = \prod_{i=1}^{N} P_i

RooFit

• pdfからNLLをつくる

RooAbsReal* nll = w.pdf("model")->createNLL(data);

最小化する

作ったNLLを最小化するには次のアルゴリズムがよく使用される

• migrad
• hesse
• minos
  前者２つは最小となる値の付近ではNLLが二次関数（パラボラ）で近似できるという仮定に基づいて計算を行う。またminosは高エネ素粒子界隈でがよく用いられる手法。minosの場合にはNLLではなく、次に定義する profile likelihood ratio（pll）がPOIとしてよく用いられる。

\lambda = \frac{L(\mu, \hat{\hat{\theta}})}{L(\hat{\mu}, \hat{\theta})}

これらの最小化コマンドは適当な順番で呼ばないとエラーが出る。

RooMinimizer min(*nll);
  {
    min.hesse(); // run Migrad before Hesse!   
  }
  {
    min.minos(); // run Migrad before Minos!   
  }

このままでは扱いづらいので、

-2\ln\lambda = -2 \left(L(\mu, \hat{\hat{\theta}}) - L(\hat{\mu}, \hat{\theta}) \right)

が使用される。こうすることでLikelihoodの引き算としてPOIを扱うことができる。

RooAbsReal* pll = nll->createProfile(*w.set("poi"))

作った pll は次のように最小化して、パラメータの推定値を算出する。

Minuit2

C++で書かれた最小化のためのライブラリ。元々のMinuitはFortranから移植されたもので、よりオブジェクト指向に寄り添って作成されたものである。

RooFit

データ点をフィットする時には $\chi^2$フィットか、$-\ln L$フィットが主流である。ROOTでは$\chi^2$、RooFitでは$-\ln L$でフィットできるという特徴がある。特に$-\ln L$フィットは統計量が小さくてもフィットできるという強みがある。

RooMinimizer

Ref

• https://root.cern.ch/download/minuit.pdf
• https://twiki.ppe.gla.ac.uk/bin/view/ATLAS/HiggsAnalysisAtATLASUsingRooStats
• https://root.cern.ch/root/html/tutorials/roofit/index.html
  • とりあえず中途半端な解説スライド見るくらいなら、tutorialコード動かすのが一番いい

RooWorkspace

ワークスペース「作業場」の名前のとおり、一旦このRooWorkspaceを作成できればあとは色々と作業を進めることができます。作成した（どうやって作成するかは後述）RooWorkspaceは通常のオブジェクトと同様Getで取ってこれます。

input_file = R.TFile("input.root") # RooWorkspace がセーブされているファイルを開く
ws = input_file.Get("workspace")   # 「workspace」という名前の RooWorkspace を取得 （もちろん適宜名前は変更）

このRooWorkspaceには、様々なオブジェクトを紐付けて保存しておくことができ、その中で最も使われるものの一つは尤度関数（Likelihood）であると思います。以下ではインプットとする分布から尤度関数を計算して、RooWorkspaceの形式に落としてくれるHistoFactoryと呼ばれるパッケージを説明します。

HistFactory

インプットとなるヒストグラムからRooWorkspaceを作成してくれる。簡単にインプットのファイルとして全１０ビンの次の様なヒストグラムを用意する。赤が信号、青が背景事象、黒点がデータ点であるような分布から尤度関数を計算→統計処理、予想感度を計算するというお決まりの流れ。

まずこのインプットファイルからRooWorkspaceを作るには次のコマンドを叩く（このコマンドはCERN ROOTライブラリ付属のもの。入っていなければ再度コンパイルしてinstallのこと）。

hist2workspace -standard_form /path/to/driver.xml

適宜XMLファイルを書くことで、どういう尤度関数を計算するかを判断してくれる。hist2workspaceが作ってくれるRooWorkspaceに保存されている内容は

- variables
- p.d.f.s
- functions
- datasets
- embedded datasets (in pdf and functions)
- parameter snapshots
- namesets
- generic objects

である。

variables

各確率密度関数の変数（「...」で冗長なoutputは省略しています）。Profile Likelihood fit をするのであれば様々なNuisance Parameterもこの欄に表示される。

variables
---------
(
Lumi,
SigXsecOverSM,
alpha_Luminosity,
binWidth_obs_x_Region_0,
binWidth_obs_x_Region_1,
channelCat,
gamma_stat_Region_bin_0,
gamma_stat_Region_bin_1,
...
nom_alpha_Luminosity,
nom_gamma_stat_Region_bin_0,
nom_gamma_stat_Region_bin_1,
...
nominalLumi,
obs_x_Region,
weightVar)

p.d.f.s

確率密度関数が定義されている。

p.d.f.s
-------
RooRealSumPdf::Region_model[ binWidth_obs_x_Region_0 * L_x_signal_Region_overallSyst_x_Exp + binWidth_obs_x_Region_1 * L_x_background_Region_overallSyst_x_StatUncert ] = 109 
RooGaussian::alpha_LuminosityConstraint[ x=alpha_Luminosity mean=nom_alpha_Luminosity sigma=1 ] = 1 
RooPoisson::gamma_stat_Region_bin_0_constraint[ x=nom_gamma_stat_Region_bin_0 mean=gamma_stat_Region_bin_0_poisMean ] = 0.0126146
RooPoisson::gamma_stat_Region_bin_1_constraint[ x=nom_gamma_stat_Region_bin_1 mean=gamma_stat_Region_bin_1_poisMean ] = 0.0132968
RooPoisson::gamma_stat_Region_bin_2_constraint[ x=nom_gamma_stat_Region_bin_2 mean=gamma_stat_Region_bin_2_poisMean ] = 0.0141033
RooPoisson::gamma_stat_Region_bin_3_constraint[ x=nom_gamma_stat_Region_bin_3 mean=gamma_stat_Region_bin_3_poisMean ] = 0.0150768
RooPoisson::gamma_stat_Region_bin_4_constraint[ x=nom_gamma_stat_Region_bin_4 mean=gamma_stat_Region_bin_4_poisMean ] = 0.0162845
RooPoisson::gamma_stat_Region_bin_5_constraint[ x=nom_gamma_stat_Region_bin_5 mean=gamma_stat_Region_bin_5_poisMean ] = 0.0178383
RooPoisson::gamma_stat_Region_bin_6_constraint[ x=nom_gamma_stat_Region_bin_6 mean=gamma_stat_Region_bin_6_poisMean ] = 0.019943
RooPoisson::gamma_stat_Region_bin_7_constraint[ x=nom_gamma_stat_Region_bin_7 mean=gamma_stat_Region_bin_7_poisMean ] = 0.0230265
RooPoisson::gamma_stat_Region_bin_8_constraint[ x=nom_gamma_stat_Region_bin_8 mean=gamma_stat_Region_bin_8_poisMean ] = 0.0281977
RooPoisson::gamma_stat_Region_bin_9_constraint[ x=nom_gamma_stat_Region_bin_9 mean=gamma_stat_Region_bin_9_poisMean ] = 0.039861
RooGaussian::lumiConstraint[ x=Lumi mean=nominalLumi sigma=0.0001 ] = 1 
RooProdPdf::model_Region[ lumiConstraint * alpha_LuminosityConstraint * gamma_stat_Region_bin_0_constraint * gamma_stat_Region_bin_1_constraint * gamma_stat_Region_bin_2_constraint * gamma_stat_Region_bin_3_constraint * gamma_stat_Region_bin_4_constraint * gamma_stat_Region_bin_5_constraint * gamma_stat_Region_bin_6_constraint * gamma_stat_Region_bin_7_constraint * gamma_stat_Region_bin_8_constraint * gamma_stat_Region_bin_9_constraint * Region_model(obs_x_Region) ] = 5.82889e-16
RooSimultaneous::simPdf[ indexCat=channelCat Region=model_Region ] = 5.82889e-16

functions

各ビンでのイベント数

functions
--------
RooProduct::L_x_background_Region_overallSyst_x_StatUncert[ 1 * background_Region_overallSyst_x_StatUncert ] = 100 
RooProduct::L_x_signal_Region_overallSyst_x_Exp[ 1 * signal_Region_overallSyst_x_Exp ] = 9 
RooHistFunc::background_Region_nominal[ depList=(obs_x_Region) depList=(obs_x_Region) ] = 100 
RooProduct::background_Region_overallSyst_x_Exp[ background_Region_nominal * background_Region_epsilon ] = 100 
RooProduct::background_Region_overallSyst_x_StatUncert[ mc_stat_Region * background_Region_overallSyst_x_Exp ] = 100 
RooProduct::gamma_stat_Region_bin_0_poisMean[ gamma_stat_Region_bin_0 * gamma_stat_Region_bin_0_tau ] = 1000
RooProduct::gamma_stat_Region_bin_1_poisMean[ gamma_stat_Region_bin_1 * gamma_stat_Region_bin_1_tau ] = 900 
RooProduct::gamma_stat_Region_bin_2_poisMean[ gamma_stat_Region_bin_2 * gamma_stat_Region_bin_2_tau ] = 800 
RooProduct::gamma_stat_Region_bin_3_poisMean[ gamma_stat_Region_bin_3 * gamma_stat_Region_bin_3_tau ] = 700 
RooProduct::gamma_stat_Region_bin_4_poisMean[ gamma_stat_Region_bin_4 * gamma_stat_Region_bin_4_tau ] = 600 
RooProduct::gamma_stat_Region_bin_5_poisMean[ gamma_stat_Region_bin_5 * gamma_stat_Region_bin_5_tau ] = 500 
RooProduct::gamma_stat_Region_bin_6_poisMean[ gamma_stat_Region_bin_6 * gamma_stat_Region_bin_6_tau ] = 400 
RooProduct::gamma_stat_Region_bin_7_poisMean[ gamma_stat_Region_bin_7 * gamma_stat_Region_bin_7_tau ] = 300 
RooProduct::gamma_stat_Region_bin_8_poisMean[ gamma_stat_Region_bin_8 * gamma_stat_Region_bin_8_tau ] = 200 
RooProduct::gamma_stat_Region_bin_9_poisMean[ gamma_stat_Region_bin_9 * gamma_stat_Region_bin_9_tau ] = 100 
ParamHistFunc::mc_stat_Region[ ] = 1 
RooStats::HistFactory::FlexibleInterpVar::signal_Region_epsilon[ paramList=(alpha_Luminosity) ] = 1 
RooHistFunc::signal_Region_nominal[ depList=(obs_x_Region) depList=(obs_x_Region) ] = 9 
RooProduct::signal_Region_overallNorm_x_sigma_epsilon[ SigXsecOverSM * signal_Region_epsilon ] = 1 
RooProduct::signal_Region_overallSyst_x_Exp[ signal_Region_nominal * signal_Region_overallNorm_x_sigma_epsilon ] = 9 

datasets

datasets
--------
RooDataSet::asimovData(obs_x_Region,weightVar,channelCat)
RooDataSet::obsData(channelCat,obs_x_Region)

embedded datasets (in pdfs and functions)
-----------------------------------------
RooDataHist::signal_RegionnominalDHist(obs_x_Region)
RooDataHist::background_RegionnominalDHist(obs_x_Region)

parameter snapshots

とある状態のパラメーターの値を保存しておくもの。尤度関数を計算していくときに、「このパラメーターはこの値に固定する」といった操作を行っていくので、初期状態であったりする値を保存しておくと便利です。

parameter snapshots
-------------------
NominalParamValues = (nom_alpha_Luminosity=0[C],nom_gamma_stat_Region_bin_0=1000[C],nom_gamma_stat_Region_bin_1=900[C],nom_gamma_stat_Region_bin_2=800[C],nom_gamma_stat_Region_bin_3=700[C],nom_gamma_stat_Region_bin_4=600[C],nom_gamma_stat_Region_bin_5=500[C],nom_gamma_stat_Region_bin_6=400[C],nom_gamma_stat_Region_bin_7=300[C],nom_gamma_stat_Region_bin_8=200[C],nom_gamma_stat_Region_bin_9=100[C],weightVar=0,obs_x_Region=9.5,Lumi=1[C],nominalLumi=1[C],alpha_Luminosity=0,gamma_stat_Region_bin_0=1,gamma_stat_Region_bin_1=1,gamma_stat_Region_bin_2=1,gamma_stat_Region_bin_3=1,gamma_stat_Region_bin_4=1,gamma_stat_Region_bin_5=1,gamma_stat_Region_bin_6=1,gamma_stat_Region_bin_7=1,gamma_stat_Region_bin_8=1,gamma_stat_Region_bin_9=1,SigXsecOverSM=1,binWidth_obs_x_Region_0=1[C],binWidth_obs_x_Region_1=1[C])

named sets

named sets
----------
ModelConfig_GlobalObservables:(nom_alpha_Luminosity,nom_gamma_stat_Region_bin_0,nom_gamma_stat_Region_bin_1,nom_gamma_stat_Region_bin_2,nom_gamma_stat_Region_bin_3,nom_gamma_stat_Region_bin_4,nom_gamma_stat_Region_bin_5,nom_gamma_stat_Region_bin_6,nom_gamma_stat_Region_bin_7,nom_gamma_stat_Region_bin_8,nom_gamma_stat_Region_bin_9)
ModelConfig_NuisParams:(alpha_Luminosity,gamma_stat_Region_bin_0,gamma_stat_Region_bin_1,gamma_stat_Region_bin_2,gamma_stat_Region_bin_3,gamma_stat_Region_bin_4,gamma_stat_Region_bin_5,gamma_stat_Region_bin_6,gamma_stat_Region_bin_7,gamma_stat_Region_bin_8,gamma_stat_Region_bin_9)
ModelConfig_Observables:(obs_x_Region,weightVar,channelCat)
ModelConfig_POI:(SigXsecOverSM)
globalObservables:(nom_alpha_Luminosity,nom_gamma_stat_Region_bin_0,nom_gamma_stat_Region_bin_1,nom_gamma_stat_Region_bin_2,nom_gamma_stat_Region_bin_3,nom_gamma_stat_Region_bin_4,nom_gamma_stat_Region_bin_5,nom_gamma_stat_Region_bin_6,nom_gamma_stat_Region_bin_7,nom_gamma_stat_Region_bin_8,nom_gamma_stat_Region_bin_9)
observables:(obs_x_Region,weightVar,channelCat)

generic objects

generic objects
---------------
RooStats::ModelConfig::ModelConfig