PCLevel
A.3 ROOT を用いたデータ解析
(付録Fを予め読む事)
ATLAS検出器で収集した実際のデータから2つのµ粒子候補を含む事象を選んで,Ntuple 形式で記録したデ
ータサンプルを用意した.全事象数は86,810である.Ntuple は事象毎に解析に必要な情報(具体的には次の
「NTUPLEに含まれる物理量」を参照)を記憶しているので,ある情報に要求を加えた(例えばpTの大きさに対
するカットを加える)場合に他の物理量の分布を表示するような処理ができ,カット値を変えて分布の変化を見
るような場合に便利である.カットや物理量の定義が単純な場合は,ROOTによる対話形式やGUIによる方式 が便利であるが,複雑な解析をするにはC++によるコードによる方が扱いやすい. wgetコマンドで以下のサン プルをダウンロードし,両方の方式で分布をみよう.
hep-www.px.tsukuba.ac.jp/~exp3/remoteClass/atlas/mumu.root
(3-1) NTUPLEに含まれる物理量:ROOT対話形式
Ntupleの構成をみよう.以下のコマンドに従い,物理量のリストを表示する.作業中にいくつかWarningが
表示されるが無視して進める.
$ root .... rootを起動する.
root[ ] TFile f("mumu.root") ... ROOT fileをオブジェクトfとして開く. root[ ] f.ls() ... ROOT fileの構造をリストする.
root[ ] TBrowser a ... 構造がGUIで表示され,ROOT Filesの下をたどればグラフ描画もできる. 表示されるグラフの横軸の適当な領域をマウスで選択するとその領域だけが拡大でき,グラフの内側,軸の上,グ ラフの外側などで右クリックすると表示方式の変更もできる.グラフの外側を右クリックして縦軸を対数表示に してみよ.GUIは右上の× で閉じる.
完成したヒストグラムの描画だけならばTBrowserで足りるが,2次元やCUTを加えた分布の表示などにはよ り多機能なTreeViewerが有効である. あるいはNtupleのオブジェクトZMASSを指定して,直接カットを加え て描画することもできる.
root[ ] TTree *t =(TTree*)f.Get("ZMASS") ... TTree形式でオブジェクト”ZMASS”をロード する.
これによりtを用いてNtupleの変数にアクセスできるようになった.
root[ ] t->StartViewer() ... TreeViewerというGUIでNtuple変数を表示する
TreeViewerの左のパネル(Current Folder)にはフォルダーが,右のパネル(Current Tree)にその中の変数 名がリストされる. 緑の葉っぱマークをクリックすると,その変数の分布が描画される. E()-emptyをクリック して適当なExpression(例えば expression=Q1*Q2<0, alias=OppCharge)を定義し,はさみマークにドロップ してから葉っぱマークを再度クリックすると,カットのかかった分布が得られる.分布や事象数に変化はあっただ ろうか.別のカットを確かめたい場合は新たにE()-emptyに条件を入れ,それをはさみマークにドラッグしなお せば良い.
2変数の相関を見る場合は,X: Y: に見たい変数をドロップしパネル左下のDrawアイコンをクリックする. また,TreeViewerのScan boxをクリックすれば,X: Y: Z:に指定したNtuple変数の値をCUTを加えた条件 でも表示できる.
TreeViewerに依らず,直接ROOT画面にDrawメソッドを入力しても良い. Draw の引数は,変数名,CUT,オ プション,イベント数,最初のイベント番号である. この方式に従って,既に定義したtオブジェクトに対して,
root[ ] t->Draw("Px1","abs(Px1)<100")
とタイプすれば,絶対値で100 GeV以下のPx1分布が得られ
root[ ] t->Draw("sqrt(Px1*Px1+Py1*Py1)","abs(Pz1)<100") root[ ] t->Draw("Px1:Py1","abs(Px1)+abs(Py1)<100")
root[ ] t->Draw("Px1:Py1","abs(Px1)+abs(Py1)<100","LEGO2")
では定義した変数や変数の相関も描画できる. “LEGO2”オプションでは,2次元プロットを色つきのLEGO表 示する.
次表に変数名とその定義を示す.ここで扱う2つのミュー粒子候補は,pTが大きい方を候補1としている.
変数名(Variable) 物理量の定義(エネルギー:GeV,長さ:cm,角度:rad)
Zvert 衝突点のZ座標(ビーム軸方向成分)
Miset 消失横向きエネルギー(missingET)の大きさ
Misph 消失横向きエネルギーのϕ方向
Px1 候補µ粒子1の運動量X成分
Py1 候補µ粒子1の運動量Y 成分
Pz1 候補µ粒子1の運動量Z成分
Q1 候補µ粒子1の電荷
Iso1 候補µ粒子1のまわり(∆R <0.2)の飛跡検出器で測定したpT和 IsoCal1 候補µ粒子1のまわり(∆R <0.3)カロリメータエネルギーET和
Px2 候補µ粒子2の運動量X成分
Py2 候補µ粒子2の運動量Y 成分
Pz2 候補µ粒子2の運動量Z成分
Q2 候補µ粒子2の電荷
Iso2 候補µ粒子2のまわり(∆R <0.2)の飛跡検出器で測定したpT和 IsoCal2 候補µ粒子2のまわり(∆R <0.3)カロリメータエネルギーET和
(3-2) Z 粒子選択のカットを決定する
Z粒子の質量を測定する場合,バックグランドが多いサンプルでは信頼度の高い測定は難しいが用意したサン プルではすでにある程度のカットがかけられているので,Z 粒子信号はカットなしでも確定しやすい. 以下には,
一般的にどのように信号を明確にするかの指針を与える.またより質量の軽い粒子を観測したい場合はバックグ ラウンドが増えるので以下の一般的な方針を参考にしてもらいたい.
1. 物理過程から分布が容易に推測できる場合は,分布を直接見てカット値を決定でき,信号を落とさないで バックグラウンドを低減できる.例えば,µ粒子は物質中でほぼ一定のエネルギーを損失するのでIsoCal は大きな値にならないが,バックグランド粒子(π,p,Kなど)は強い相互作用で反応し,大きな信号とな る. 本来はシミュレーションで最適値を決めるが,分布の不連続牲から判断できる場合もある. 分布から外 れるものはZ粒子信号領域の外に拡がる「裾」(分布のtailと呼ぶ)の成分となるので,カットの設定が適
切かの判断ができる.
カットを与えないと,ROOTは,含まれるデータをすべて表示するようヒストグラムの範囲を決めてしま う.検出器のノイズやバックグラウンドなどにより意味の無いデータがあると,本当に見たいµ粒子の 信号領域が狭められて,1 ∼2 binに押し込められてしまうこともある.例えば,陽子陽子衝突は8 TeV であるためこれを超えた運動量を持つµ粒子は発生しない.Z粒子の質量を考えると崩壊して発生する µ粒子の運動量成分はせいぜい100 GeVである.極端に大きな運動量のものは宇宙線のバックグランド と考えられる. そのため,既述のカットを加えた条件で分布を見るようにしよう. 信号領域の目安として
「A.2ATLASでのµ粒子の同定」で示されているカット値を参考にする.
2. 信号領域があいまいな場合は,他のカットを加えた上での分布をみて決定する.他のカットによりバック グランドが低減され,信号領域の分布が明らかになることが期待できる.
3. それでも判別しない場合はZの分布のピーク部分のイベントが減らないようする. Zの分布は次節で扱う のでカットは緩めな値に仮設定しておく.
(3-3) Z粒子の質量
カットを加えた上でZ粒子の不変質量を組むためには,ROOTのViewerを用いても可能ではあるが,ここで はC++でプログラムを書いて実行しよう. C++でNtupleデータを読み出すには,ROOTデータファイルの 構造を受け渡さなくてはいけない. ROOTはそのために必要なヘッダーファイルとイベントループをするだけで 何もしないC++ファイルを生成してくれる .これらのファイル(Zmass.C, Zmass.h)にヒストグラムを定義 する部分を書き加えたファイルmyZmass.Cを用意したので,それらを自分の作業領域にダウンロードしよう.
$ wget http://hep-www.px.tsukuba.ac.jp/~exp3/remoteClass/atlas/myZmass.h .
$ wget http://hep-www.px.tsukuba.ac.jp/~exp3/remoteClass/atlas/myZmass.C .
プログラムを実行するために,ROOTを再起動して root[ ] .L myZmass.C... myZmass.Cをロードする root[ ] myZmass Z ... myZmassからオブジェクトZを生成 root[ ] Z.Loop()... myZmassのLoopメソッドを実行
Terminal上で$ emacs myZmass.C &としてEditorを起動して,myZmass.Cを編集しよう. Zの不変質量を計算するために,関数文 massをプログラムの最初に用意している. 必要な変数を渡 して,正しく不変質量を返すように編集せよ. 今は仮の値として90を返すようにしている.
いくつかの変数に対しては仮のカットの値が設定してある. 最適と思われるカット値に設定せよ. ま た,これ以外の変数に対する要求があるならば,適切にプログラムを変更せよ.
全体のCUTはいくつかの定義済のBoolean変数 の論理式で定義してある. 変数を加えた場合は,
$myZmass.C中の以下のCUTの定義式の部分に,それらを正しくに反映させなくてはならない.
節(3-1)で行ったようにTTreeオブジェクトtを生成した後にt->MakeClass(”Zmass”)とすると,Zmass.C, Zmass.hができる.
true(真)かfalse(偽)の二値を取る論理変数の型.if(CUT)の条件文はCUTがtrueの場合に実行される.尚,CでのAとBの論 理積(AND)はA && B,AとBの論理和(OR)はA||B,Aの否定(NOT)は!Aである.数値の比較には,等号(=)は==,そ の否定(�)は!=,大小関係(<,>, �, �)には,それぞれ<,>,<=,>=を使う.
bool CUT = Cut1 && Cut2 && Cut3;
用意したヒストグラム以外の分布を見たい場合は,相当するヒストグラムを用意して描画するように 変更せよ.ヒストグラムはCUTを与える前後で分布を比べられるようにしている.
hC->Draw() .... CUTを加えたヒストグラムを描画 h->Draw(“same”) .... CUT前のヒストグラムを重ね描き
Zの分布が得られたら,ガウス関数をフィットして中心値を求めてみよう. カットを加えたZ分布をフィット するために,ヒストグラムオブジェクトh_ZmassFを用意してあるので,ROOTを起動してmyZmassのループメ ソッド実行後にガウスフィットする.
root[ ] Z.Loop()... ここまでは既述 root[ ] h_ZmassF->Fit("gaus","","E1") (色の設定や再描画は適度にやる)
Z粒子の不変質量分布はガウス関数では再現できないことが分かるであろうか. 合わない領域を無視してピー クの周辺のみをフィットする時は,
root[ ] h_ZmassF->Fit("gaus","","E1",lower,higher)
とする. ここでlowerなどはフィットする領域を限定するための横軸の値である.
Z 粒子の質量はガウス分布の中央値としてフィットの結果から求められる.rootの画面にはフィット結果 が示される.フィット結果をグラフにも表示するためにはグラフ右上の統計数を示すボックス内で右clickし SetOptFitwindowを表示させ,値とし1111を入力する.
中央値のタイプAの不確かさは,ガウス分布の標準偏差σとピークに含まれる事象数N より,
∆M =σ/√ N
で与えられる.フィットの結果が与える標準不確かさと比較してみよう.
きついカットを適用するとイベント数が少なくなり統計のばらつきが増すので,必ずしも良い測定結果になら ない. 例えば,2つのミュー粒子の片方はきつくしても,もう片方にはゆるめのカットを適用することも可能であ る. カット値を調整してみよう.
バックグラウンドの少ない場合(Zの場合)はN としてピークに含まれる事象数N を取れば良いだろう. しか し,例えばJ/ψの場合,ピークにはバックグラウンドの数も含まれる. この場合,タイプA(統計)の不確かさ はどうなるであろうか考察せよ. 最終的には,タイプAの不確かさに加えて,タイプBの不確かさも考慮しなく てはならないが,その評価には検出器の細部に至る理解が不可欠であり,授業では省略する.
Z粒子はCERNのLEP電子・陽電子衝突器で精密測定されている.2010年版のParticle Dataによると質 量は91.1876±0.0021GeV/c2 である.測定結果をこの値と比較せよ.
■軽い粒子(J/ψ)の質量 ミュー粒子対に崩壊する粒子はZ以外にもある.1974年に発見されたc¯cの束縛状態 であるJ/ψ(3096.916±0.011MeV/c2)や1977年に発見されたb¯bの束縛状態であるΥ(Υ(1S)は9460.30±0.26 MeV/c2)も各々6.0%,2.5%の崩壊分岐比でミュー粒子対に崩壊する.mumu.rootのデータでJ/ψはどの様に 見えるか.