Documents for ATLAS-Japan Software Tutorial

Introduction

このページではATLAS(に入る予定)の人のためにC++とROOTあたりからプログラミング・解析のためのTipsをまとめています。
ただし可能な限りATLASに特化しない方向でいきたいと思います。なので別の実験に移っても有用なことしかやりません。

C++ tutorial

ここでは教科書に書いてあるような基本的なことはスキップします。
世の中にはたくさんの教科書がありますし、今の時代ググれば基本情報は出てくるので、以下は特に知っておいた方が良いことを書きます。
(+参加者だけはホワイトボードでオブジェクト指向の基本をおさらいします。)

したがってそこで不安になる方は

あたりで勉強することをお勧めします。

知っておくとためになる知識

オブジェクトを作る(消す) - memory allocation, ポインタ, 参照
HELP Will Buttinger - offline software tutorial 2015のスライドより拝借しています

オブジェクトをC++で作ると、それ用のメモリーがあてがわれます。 これには2種類 (stackheap)があります。

MyClass obj;   // creating instance of MyClass on the stack
int myInt;       // creating instance of int on the stack
MyClass* obj2 = new MyClass();    // creating instance of MyClass on the heap
int* myInt2 = new int();              // creating instance of int on the heap 
stackは
  • read/writeが速い
  • heapよりはサイズが小さい
  • 自動的に削除される(自分で delete する必要がない)
一方でheapは
  • stackよりは遅い
  • stackよりはサイズが大きい
  • 自分で delete などをしないといけない(ROOT上では多少やってくれますが)

一つ例を上げると

int main()
{
  int a;
  int* b = new int;     // note: no brackets means ‘using default constructor’
  {
    int c;
    int* d = new int;
  }      // At this point c is deleted. The int that d points to hasn’t though... MEMORY LEAK!!
  delete b;      // good, we’re cleaning up our heap allocation
}     // a gets deleted from the stack here.

ポインタとは一言で言えば、メモリ上のインスタンスのアドレスです。

MyClass obj;    // creating instance of MyClass on the stack
int myInt;        // creating instance of int on the stack
MyClass* obj2 = new MyClass();     // creating instance of MyClass on the heap
int* myInt2 = new int();              // creating instance of int on the heap

MyClass* obj3 = &obj;     // obj3 points at obj
int* myInt3 = &myInt;     // myInt3 points at myInt

アドレスってのはメモリではとても小さいです (64bitマシンで8 bytes (=64 bits))
なので関数に受け渡すときなどはオブジェクト自体を渡すよりもポインタを渡したほうが早くなります。

  • basicなtype (bool, int, float, ...)では差がありませんが。
  • 関数の引数でインスタンスを与えると同じ内容のオブジェクトがstack上にコピーされます。
    そしてその関数が終わるときに自動的に消されます。(i.e. オリジナルのオブジェクトにはアクセスしていません)
  • 引数にポインタを与えるとアドレスのコピーが作られます。
    コピーされたアドレスにアクセスすることになるのでオリジナルのオブジェクトにアクセスすることになります。
void myBadFunction(MyClass obj) { obj.myMethod(); }
void myGoodFunction(MyClass* obj) { obj->myMethod(); }

int main()
{
  MyClass obj;
  myGoodFunction(&obj);
}

なんかしらオブジェクトを作ったら、必ず「いつ」「どこで」消える(消す)のかを考える癖を付けましょう。
stackの場合、自動で削除されます。でもscope (中括弧 { }) の外側ではアクセスしてはいけません。

MyClass* b;     // created a pointer, but not made it point at anything
{
    MyClass c;
    b = &c;    // b now points at c.
}    // c is deleted, but we’ve dangerously kept it’s address (held in b)
b->myMethod();    // Can’t do this, c doesn’t live at the address held by b anymore

heapの場合: 必ず自分でdeleteをしないといけません。
ダメな例を示すと

void myBadFunction() {
     MyClass* obj = new MyClass();
     if(obj->someMethod()) return;
     delete obj;     // too late if someMethod() returned true...
}

最後に参照について復習しましょう。

  • MyClassタイプのインスタンスのポインタの型は MyClass* です。
  • MyClassタイプのインスタンスの参照の型は MyClass& です。
  • 参照はアドレスの代わりにオブジェクトを指すことができます。
MyClass obj
MyClass* objPtr = &obj;
MyClass& objRef = obj;

継承, Virtual 関数, キャスト

メンバー関数を一つ持っているクラスと、それを継承したクラスを作ってみます。

class MyClassA
{
 public:
  void myMethod() { std::cout << "It's MyClassA" << std::endl; }
};

class MyClassB : public MyClassA
{
 public:
  void myMethod() { std::cout << "It's MyClassB" << std::endl; }
};
これを実際に使用してみるとポインタ / 参照した型のタイプに依存した結果が示されます。
MyClassB* b = new MyClassB();
MyClassA* b_as_a = b;

b->myMethod();        // prints "It's MyClassB"
b_as_a->myMethod();   // prints "It's MyClassA"
ここで関数をvirtualにしてみましょう。
class MyClassA
{
 public:
  virtual void myMethod() { std::cout << "It's MyClassA" << std::endl; }
};

class MyClassB : public MyClassA
{
 public:
  virtual void myMethod() { std::cout << "It's MyClassB" << std::endl; }
};
するとコンパイラが呼ぶ関数は実際のインスタンスのものになります。
MyClassB* b = new MyClassB();
MyClassA* b_as_a = b;

b->myMethod();        // prints "It's MyClassB"
b_as_a->myMethod();   // prints "It's MyClassB"
このようなことをポリモーフィズム(多様性)と言います。

キャストについても学んでおきましょう。

  • C-style cast
  • static_cast
  • dynamic_cast
  • reinterpret_cast
  • const_cast
の5つありますが知っておいてほしいものはstatic_castとdynamic_castだけです。それ以外は極力使わないようにしましょう。細かいことは言いませんが危険をはらんでいるんです。

dynamic_cast

安全な基底クラスから派生クラスへのキャスト(ダウンキャスト) のために使います。
(派生クラスから基底クラスへのキャスト(アップキャスト)も特に問題なくできるが、そちらはstatic_cast推奨)

static_cast

主に派生クラスから基底クラスへのキャスト(アップキャスト)用。皆さんがやりたい(と思われる)C-styleキャストはこれに置き換えましょう。
危険なキャストをコンパイルエラーではじく事は可能です。

使い方一つだけ例を示すと

int ival;
long lval = static_cast<long>(ival);

Coding rule

言いたいことは実はここだけであったりするくらい重要なお話です。とにかく
  • 自分が読んですぐに理解できるコードを書く!
  • 他の人も使う(使って欲しい)ならば、シェアしうるコードを書く!
とにかく一定のルールに従ったコードを書くことが重要です。(例えばメンバー変数なら m_ で始めるとか)

ATLASでは推奨コーディング法というものがあって、ここに最新のがあります。

一部古い情報もありますが日本語版も坂本先生が用意してくれています。

あと、一つのソースコードのファイルが1000行とか超えるようならそれはデザインが間違えてると思いましょう。 きちんとしたオブジェクト志向にのっとって作っておけば(コードの数は多少増えても)そういうことは起こらないはずです。
(1万行を超えるコードなんて読む気も失せるでしょ?)

Exercise 1 - C++ programming

自分で一つクラスを用いたプログラムを作ってみましょう。
ここでは
  • C++コードの基本的なコンパイル方法を知る
  • 一つのプログラムのために複数のコードを用いる(ヘッダファイルを作る)
  • クラスを作ってみる
にフォーカスします。

とりあえずどこか作業する場所を作ります。

[junpei@login01 ~]$ mkdir cxxtut2015
[junpei@login01 ~]$ cd cxxtut2015

ではParticleという粒子のクラスを作ります。 このクラスは

  • 粒子の名前
  • エネルギー
  • 電荷
の3つのパラメータを保持するとします。また粒子の情報を出力できるという機能を持つことにします。 クラス定義はヘッダファイル (Particle.h) を書いていきましょう。
#ifndef Particle_h
#define Particle_h
 
#include <string>
 
class Particle
{
 public:
  Particle();
  virtual ~Particle();
 
  virtual void printInfo();

  inline void setName(std::string name) { m_name = name; }
  inline void setEnergy(double energy) { m_energy = energy; }
  inline void setCharge(const int charge) { m_charge = charge; }
 
  inline std::string getName() const { return m_name; }
  inline double getEnergy() const { return m_energy; }
  inline int getCharge() const { return m_charge; }
  
 protected:
  std::string m_name;             ///< particle name
  double m_energy;             ///< 4-momentum
  int m_charge;                           ///< charge
 };
 
#endif

Particleクラスの関数の中身は(Particle.cxx)に実装しましょう。

#include "Particle.h"

#include <iostream>

Particle::Particle()
{
  m_name = "NOT DEFINED";
  m_energy = 0.0;
  m_charge = 0;
}
 
Particle::~Particle()
{}

void Particle::printInfo()
{
  std::cout << "Particle Name = " << m_name
           << ", energy = " << m_energy
           << ", charge = " << m_charge << std::endl;
}

そしてmain()を作ってクラスを呼びます (exercise.cxx)

#include "Particle.h"

int main(void)
{
  Particle electron;

  electron.setName("electron");
  electron.setEnergy(1000.);   // 1000 MeV
  electron.setCharge(-1);

  electron.printInfo();

  return 0;
}

ではコンパイルしてみます。

[junpei@login01 cxxtut2015]$ g++ -c Particle.cxx
[junpei@login01 cxxtut2015]$ g++ exercise.cxx Particle.o -o exercise.exe

せっかくなので実行
[junpei@login01 cxxtut2015]$ ./exercise.exe

ROOT tutorial

ROOTとは

一般的に我々の実験では

  1. 何か測定したい対象 → 検出器・システムを構築する
  2. そのシステムからの情報を収集・加工する
  3. そのデータをさらに解析し、最終的な結果を得る
ROOTはそのうち2, 3のために使われる「フレームワーク」で、 具体的には物理計算、ヒストグラムやグラフの作成、フィッティング、・・・を行うことが出来るツールである。

どういうものを "プロット" するのか?

fig_02a.png atlasandcmse.png

ヒストグラムとTree (n-tuple)

ヒストグラムとTreeについて、先にざっと説明してしまいます。
ヒストグラム
ヒストグラムとは試行したら値がXだった、次はYだった、さあ(確率)分布は? というものです。つまり
  • 値=X → Xの該当するビンを+1する
ということが基本になります。

ヒストグラムに限らず、ROOTではたくさんの「もの」に名前がついています。
hist.png
これ知っていると(いずれ綺麗なプロット作るときに)めちゃくちゃ便利です。

ROOTではこのヒストグラムを作るためのクラスが用意されています。

TH1_Inh.gif

TH1 は全てのヒストグラムにとっての基底クラスです。ROOTでいう「もの」は全て TObject クラスを継承しています。
上の絵にあるクラスの説明をすると

  • TNamed: ROOTオブジェクトは名前を管理されている。そのためのbase class
  • TAtt[Line,Fill,Marker]: 飾り付け(お絵かき)のクラス
  • TH[1,2,3]*: 1次元、2次元、3次元
  • TH*[S,F,D...]: Short, Float, Doubleの頭文字。Bin contentsの最大値とかが違う

Tree (n-tuple)

実験データは(一般的に)数多くの情報を持っている。

  • エネルギー、運動量、位置、電荷、・・・などの粒子ごとの情報
  • これらは事象ごとには独立、でも1事象内は結びついている

例えば

荷電粒子の軌跡No 横運動量 [MeV] 座標 z [mm] 電荷
0 3200.43 7643 -1
1 2893.22 9834 1
2 3603.9 11232 1
3 9899.14 10232 -1
4 5674.32 8092 1
5 3432.33 7662 1
この1行の組をtuple (タプル or テュープル) と、複数からなるもの(表)をNtuple (エヌタップル or エヌテュープル)と呼ぶ。

ROOTではこれらを扱うものにTNtuple, TTreeといったクラスがあります。

  • TNtuple: 基本的に数字のみを扱う
  • TTree:数字だけでなく、色々な情報(文字列、STL container, TObject, ...)を保持できる (基本こっち)。

HELP 言葉の定義ですが表の横の項目を Branch 、一行のかたまりを Entry と言います。(なので複数行のことは Entries

Exercise 1: Interactiveに遊んでみる

手を動かしてやってみましょう。

まずはlogin.icepp.jpでのおまじないをして、その後作業場所を作ります。
[junpei@login01 ~]$ setupATLAS
[junpei@login01 ~]$ lsetup root
[junpei@login01 ~]$ mkdir ROOTtut2015
[junpei@login01 ~]$ cd ROOTtut2015

起動は
[junpei@login01 ROOTtut2015]$ root
スプラッシュが出た後、CINTというものが走るプロンプトが現れて起動しましたね?

終了するときは .q (ドットと小文字のQ) で終了できます
root [0] .q

HELP TIPS

  • 最初のスプラッシュとかがウザいという人には
    [junpei@login06 ~]$ root -l
    root [0]
    といきなり起動する。
  • グラフィックが必要ない (X飛ばしたくない) 場合は
    [junpei@login06 ~]$ root -l -b
    root [0]
    とバッチモードで起動する
  • 実行スクリプト(ROOTマクロ)を実行するのなら引数で与えられる
    [junpei@login06 ~]$ root -l test.C
    root [0]
    Processing test.C...
    root [1]
  • マクロ実行後、そのままROOTを終了して良いのなら
    [junpei@login06 ~]$ root -l -b -q test.C
    root [0]
    Processing test.C...
    [junpei@login06 ~]$
    となる
  • root上では . (ドット) で始まる文字列はコマンドとして解読される。それ以外はC++としての命令だと解釈される
  • .! はシェルエスケープします。例えば
    [junpei@login01 ~]$ ls
    2012_10_03S  setup.sh  test.C  tutorial2014
    [junpei@login01 ~]$ root -l
    root [0] .!ls
    2012_10_03S  setup.sh  test.C  tutorial2014
    他にも .x (マクロ実行)とか .L (マクロ読み込み)とかあります。

ヒストグラムを作る
それではとりあえずヒストグラムを作ってみましょう
[junpei@login01 ROOTtut2015]$ root -l
root [0] TH1F* myHist = new TH1F("myHist", "My first histogram", 10, 0., 10.);
TH1FクラスのオブジェクトmyHistを作ります。
とにかくこれをDrawしてみましょう
root [1] myHist->Draw();
新しいWindowが出て軸だけのものが表示されましたか?
emptyhist.gif
TH1Fのコンストラクタの引数が表示されたもののどこに対応しているか確認しましょう。

ヒストグラムにFillする
先程説明しましたがヒストグラムの基本は
  • 値=X → Xの該当するビンを+1する
です。それを実現するのが TH1::Fill() です。例えばX=3.4なら…
root [2] myHist->Fill(3.4);
もう一回Drawすると絵が更新されます。
root [3] myHist->Draw();
今、0~10を10ビンに分けたヒストグラムを作っているので、X=3~3.9999999...が同じビンに入ることになります。

HELP TIPS

  • Fill関数に引数の違いで何パターンかあります。例えばTH1クラスでは2個めの引数で重みをつけることができます(普段は1)
virtual Int_t Fill(Double_t x);
virtual Int_t Fill(Double_t x, Double_t w);
  • ただし、x=3を50回fillしたのと、1回をw=50でfillしたのは全然意味が違います。
    → そのビンの統計エラーを考える!
    • TH1は何回Fillしたか(entries)とか色々な情報を持っているので、意味が違うことはやらないこと
  • ではどういう時に使うのか?
    • MCの分布を作る (一般的にMCは大量のイベントを生成させてLumiでノーマライズする)
    • サンプルや場合によってWeightが変わるときなど・・・
  • 基本的に引数1個増やすとweightになる
    c.f.) TH2 (2次元ヒスト) クラスの場合
virtual Int_t Fill(Double_t x, Double_t y);
virtual Int_t Fill(Double_t x, Double_t y, Double_t w);

飾り付け
TH1クラスはTAttLineとかのクラスを継承しているので、その部分の飾り付けが出来ます。 基本的に Set[Line,Fill,Marker][Color,Style,...](...) というメンバ関数を使います。
root [*] myHist->SetLineColor(2);
root [*] myHist->SetFillColor(5);

TH1オブジェクトはTAxisオブジェクトを保持していて、こいつが軸の情報を持っています。 なので...

root [*] myHist->GetXaxis()->SetTitle(“abcde”);
root [*] myHist->GetYaxis()->SetLabelSize(0.03);
とかで飾り付けができます。

カラーコード: https://root.cern.ch/root/html/TColor

  • もしくはCanvasにあるメニューバーから View → Colors

とにかくクラスリファレンスが最強! むしろこれさえ使えれば自分で調べられます。

実際にはgoogleでROOT TH1とかやればすぐにそのページに飛べます。

統計情報の変更は 

root [*] gStyle->SetOptStat(1111211); // Mean with error&#12289;integral, under/overflow
root [*] gStyle->SetOptStat(0);   // no statistics information
この定義は
The parameter mode can be = ksiourmen (default = 000001111)
  k = 1; kurtosis printed
  k = 2; kurtosis and kurtosis error printed
  s = 1; skewness printed
  s = 2; skewness and skewness error printed
  i = 1; integral of bins printed
  i = 2; integral of bins with option "width" printed
  o = 1; number of overflows printed
  u = 1; number of underflows printed
  r = 1; rms printed
  r = 2; rms and rms error printed
  m = 1; mean value printed
  m = 2; mean and mean error values printed
  e = 1; number of entries printed
  n = 1; name of histogram is printed
となっています。

先に述べておくとFitの情報はgStyle->SetOptFit(Int_t)で行えます。

The parameter mode can be = pcev  (default = 0111)
    p = 1;  print Probability
    c = 1;  print Chisquare/Number of degress of freedom
    e = 1;  print errors (if e=1, v must be 1)
    v = 1;  print name/values of parameters

TLegend (Legend, ヒストグラムの説明を表示する)

root [*] TLegend* tl = new TLegend(0.52, 0.72, 0.72, 0.82);
root [*] tl->AddEntry(myHist, “tutorial example”, “F”);
root [*] tl->SetFillColor(0);
root [*] tl->Draw();

TLine (線を書く)

root [*] TLine* line = new TLine(4.0, 0.0, 6.0, 0.5);
root [*] line->SetLineColor(4);
root [*] line->Draw();

保存/読み込み
最後にせっかくやったものは保存しましょう。

絵として保存する

root [*] c1->Print(“test.pdf”);
c1 というのは TCanvas* オブジェクト (Draw() したときに自動で作られたものです)

ヒストグラムなど、オブジェクトを保存する ROOTファイルというのはその気になればTObjectクラスを継承したものなんでも保存できる

root [*] TFile* file = new TFile(“test.root”, “RECREATE”);
root [*] file->cd();
root [*] myHist->Write();
root [*] file->Close();

開く、そしてmyHistをもう一回Drawする

root [*] TFile* file = TFile::Open(“test.root”, “READ”);
root [*] TH1F* myHist2 = dynamic_cast<TH1F*>(file->Get(“myHist”));
root [*] myHist2->Draw();

マクロを使う

今までやったことを毎回やるのは面倒ですね。例えば失敗したり、微妙に変更したかったら一からやり直しになってしまいます。
というわけで一連の命令をあらかじめ別ファイルで記述するようにしましょう。

(基本的に)C++の文法に準拠するファイルを用意します。
もし一つの関数内で実行できる場合

{
#include <iostream>
  TH1F* myHist = new TH1F(“myHist”, “test histogram”, 10, 0, 10);
  for(int i=0; i<10; i++) {
    myHist->Fill(i%3);
  }
  myHist->Draw();
  std::cout << “done” << std::endl;
}
みたいな感じで書けます。 しかし現実的なことを考えると、推奨は 下のようなもので、これを test.C と保存することにします。
#include <iostream>

void HistAttribute(TH1* hist)
{
  hist->SetLineColor(2);
}

void test()
{
  TH1F* myHist = new TH1F(“myHist”, “test;x;y”, 10, 0, 10);
  TH1F* myHist2 = new TH1F(“myHist2”, “test 2”, 10, 0, 10);
  HistAttribute(myHist2);
  for(int i=0; i<10; i++) {
    myHist->Fill(i%3);
    myHist2->Fill(i);
  }
  myHist->Draw();
  myHist2->Draw(“SAME”);
  std::cout << “done” << std::endl;
}
のように関数が使えるようにしておくべきです、
ここで注意しておいて欲しいことは ファイル名の関数が必要(実行時にそいつが呼ばれる。 main() ではない) ということです。

これをROOTで実行するには
[junpei@login02 ROOTtut2015]$ root -l test.C
となります。 test.C なので test() が呼ばれます。

この test.C の中にも一部重要なポイントが入っています。口頭で説明します。

Exercise 2: Histogram and Fitting

一つ練習問題をやってみましょう。上では直接教えていないことも少しだけ入ってます。
  • 自分で平均値0, sigma=1のガウス分布になる乱数を50万事象発生させる
  • それをヒストグラムに詰める
  • これをガウス分布で-2~2の範囲をFitする
  • 自分で好きなサイズのCanvasを作る
  • 軸やLegendをつける
  • ヒストグラムをROOTファイルで保存する

これを行うマクロをexercise2.Cとしましょう。 各ステップ実装するたびに実行してみてどうなるか見ると勉強になると思います。

0. マクロファイルを用意してヒストグラムを作る
まずは枠組みを用意しましょう。

void exercise2()
{
  TH1F* hist = new TH1F("hist", ";X;Events", 100, -10, 10);

}

1. ヒストグラムを作る、ガウス分布になる乱数を発生させてFillして、Drawする

今回ガウス分布の乱数を発生させるには、 TRandom3 というクラスを利用します。
(自力でやりたい方はボックス=ミュラー法とか勉強して下さい)

void exercise2()
{
  TH1F* hist = new TH1F("hist", ";X;Events", 100, -10, 10);

  TRandom3 rndm(0);
  for(int i=0; i<500000; i++) {
    float x = rndm.Gaus(0, 1);
    hist->Fill(x);
  }
  hist->Draw();
}

2. キャンバスを自分で作る、ヒストグラムの飾り付け、Fitする
Fitに関してはROOT User's guide - Fitting Histogramsを読んで、ちゃんと勉強することをおすすめします。
(今日は時間ないのでカバーしません)

void exercise2()
{
  gStyle->SetOptStat(0);

  TCanvas* c1 = new TCanvas("c1", "", 800, 540);
  c1->Draw();
  c1->SetTopMargin(0.03);
  c1->SetBottomMargin(0.14);
  c1->SetRightMargin(0.03);
  c1->SetLeftMargin(0.14);
  c1->SetTicks();

  TH1F* hist = new TH1F("hist", ";X;Events", 100, -10, 10);
  hist->GetXaxis()->SetTitleSize(0.06);
  hist->GetYaxis()->SetTitleSize(0.06);

  TRandom3 rndm(0);
  for(int i=0; i<500000; i++) {
    float x = rndm.Gaus(0, 1);
    hist->Fill(x);
  }
  hist->Draw();

  TF1* gaus = new TF1("f1", "gaus");
  hist->Fit("f1", "", "", -2.0, 2.0);
}

3. Legendつける
void exercise2()
{
  gStyle->SetOptStat(0);

  TCanvas* c1 = new TCanvas("c1", "", 800, 540);
  c1->Draw();
  c1->SetTopMargin(0.03);
  c1->SetBottomMargin(0.14);
  c1->SetRightMargin(0.03);
  c1->SetLeftMargin(0.14);
  c1->SetTicks();

  TH1F* hist = new TH1F("hist", ";X;Events", 100, -10, 10);
  hist->GetXaxis()->SetTitleSize(0.06);
  hist->GetYaxis()->SetTitleSize(0.06);

  TRandom3 rndm(0);
  for(int i=0; i<500000; i++) {
    float x = rndm.Gaus(0, 1);
    hist->Fill(x);
  }
  hist->Draw();

  TF1* gaus = new TF1("f1", "gaus");
  hist->Fit("f1", "", "", -2.0, 2.0);

  TLegend* legend = new TLegend(0.7, 0.7, 0.95, 0.95);
  legend->AddEntry(hist, "Gaus(0, 1)", "F");
  legend->AddEntry(gaus, "Fit", "L");
  legend->SetBorderSize(0);
  legend->SetFillColor(0);
  legend->Draw();
}

4. 保存する
void exercise2()
{
  gStyle->SetOptStat(0);

  TCanvas* c1 = new TCanvas("c1", "", 800, 540);
  c1->Draw();
  c1->SetTopMargin(0.03);
  c1->SetBottomMargin(0.14);
  c1->SetRightMargin(0.03);
  c1->SetLeftMargin(0.14);
  c1->SetTicks();

  TH1F* hist = new TH1F("hist", ";X;Events", 100, -10, 10);
  hist->GetXaxis()->SetTitleSize(0.06);
  hist->GetYaxis()->SetTitleSize(0.06);

  TRandom3 rndm(0);
  for(int i=0; i<500000; i++) {
    float x = rndm.Gaus(0, 1);
    hist->Fill(x);
  }
  hist->Draw();

  TF1* gaus = new TF1("f1", "gaus");
  hist->Fit("f1", "", "", -2.0, 2.0);

  TLegend* legend = new TLegend(0.7, 0.7, 0.95, 0.95);
  legend->AddEntry(hist, "Gaus(0, 1)", "F");
  legend->AddEntry(gaus, "Fit", "L");
  legend->SetBorderSize(0);
  legend->SetFillColor(0);
  legend->Draw();

  TFile* file = new TFile("output.root", "RECREATE");
  file->cd();
  hist->Write();
  file->Close();
}

Exercise 3: Create a basic x-y plot

login.icepp.jp:/home/junpei/ROOTtut2015/Exercise3/data.txt にはxの値、yの値、yの誤差という順で並んだtext fileがあります。
このデータを使って以下のようなプロットを自分で作ってみましょう。
exercise2.gif
このExerciseでわかってほしいポイント(テクニック)は
  • 誤差付きのグラフを描く
  • 点のスタイル(Marker)を変更する
  • X, Y軸にタイトルを付ける、大きくする。
  • X軸は対数表示にする
  • X軸Y軸の表示する範囲を自由に設定する
  • キャンバスの中でプロットを表示する場所を調整する
  • エラーバーの先っちょをの線を無くす
などです。
まずはグラフオブジェクトを作る、テキストファイルを読み込む、表示してみる。

誤差付きのグラフは TGraphErrors というクラスを利用します(クラスリファレンス)。

使い方は色々あるのですが、使いやすい(と僕が思う)ものに TGraph::SetPoint(Int_t i, Double_t x, Double_t y)TGraphErrors::SetPointError(Int_t i, Double_t x, Double_t y) があります。

{
  std::ifstream fin("data.txt");
  double x, y, yerr;

  TGraphErrors* g1 = new TGraphErrors();

  int myindex = 0;
  while(fin >> x >> y >> yerr) {
    g1->SetPoint(myindex, x, y);
    g1->SetPointError(myindex, 0, yerr);
    myindex++;
  }
  fin.close();

  g1->Draw("APE");
}
exercise2-1.gif

飾り付けしていく
ベースは出来たので、あとは一つずつ飾り付けしていきます。
  • Markerのスタイルや色を変更する
    TAttMarker::SetMarkerStyle() とか
  • Canvasを独自に用意してTick, Grid, 余白のスペースを指定する
    TCanvas を使ってキャンバスを先に作っておきましょう。そうするとそいつにGrid描け、Tick付け加えろ、とか色々命令できるようになります。
  • X軸を対数表示にする
    TPad::SetLogx() を使う。
  • 好きな範囲でグラフを表示する
    TGraphはDrawオプションに "A" を付けることで自動的に軸が描かれます。
    これは多少テクニックが必要で、まず空の(2次元)ヒストグラムを書きます。それを描いて、その中に "A" オプション無しでグラフを書くということで可能となります。
  • エラーバーの最後のちょこっとした線を無くす
    TStyle::SetEndErrorSize(Float_t) という関数があります。実際には gStyle->SetEndErrorSize(0); のように使います。
全部実装すると…

{
  gStyle->SetOptStat(0);
  gStyle->SetEndErrorSize(0);

  ifstream fin("data.txt");
  double x, y, yerr;

  TGraphErrors* g1 = new TGraphErrors();

  int index = 0;
  while(fin >> x >> y >> yerr) {
    g1->SetPoint(index, x, y);
    g1->SetPointError(index, 0, yerr);
    index++;
  }
  fin.close();

  g1->SetMarkerStyle(21);
  g1->SetMarkerColor(4);
  g1->SetLineColor(4);

  TCanvas* c1 = new TCanvas("c1", "");
  c1->Draw();
  c1->SetLogx();
  c1->SetGrid();
  c1->SetTicks();
  c1->SetRightMargin(0.04);
  c1->SetTopMargin(0.03);
  c1->SetBottomMargin(0.11);
  c1->SetLeftMargin(0.11);

  TH2D* hAxis = new TH2D("hAxis", ";setting value;measured value",
                          1, 0.01, 10000, 1, 10, 30);
  hAxis->Draw("AXIS");
  hAxis->GetXaxis()->SetTitleSize(0.05);
  hAxis->GetYaxis()->SetTitleSize(0.05);
  g1->Draw("PE");
}

Treeを読む

/home/junpei/ROOTtut2015/myFirstTree.root に上でTreeについて説明した表と同じ内容のものを用意しました。
自分のところにコピーしてきたら、色々みてみましょう。
[junpei@login02 ROOTtut2015]$ root -l myFirstTree.root 
root [0] 
Attaching file myFirstTree.root as _file0...
(class TFile *) 0x3c8b4f0
root [1] .ls                                                                  → 何が入っているのか確認
TFile**         myFirstTree.root
 TFile*         myFirstTree.root
  KEY: TTree    myFirstTree;1   My First Tree for Tutorial                    → myFirstTreeというTTreeオブジェクトが入っている
root [2] myFirstTree->Print();                                                → Treeの情報を見る。定義されている変数が一覧される
******************************************************************************
*Tree    :myFirstTree: My First Tree for Tutorial                             *
*Entries :        6 : Total =            2256 bytes  File  Size =        943 *
*        :          : Tree compression factor =   1.02                       *
******************************************************************************
*Br    0 :track_pt  : track_pt/D                                             *
*Entries :        6 : Total  Size=        624 bytes  File Size  =        130 *
*Baskets :        1 : Basket Size=      32000 bytes  Compression=   1.00     *
*............................................................................*
*Br    1 :track_z   : track_z/D                                              *
*Entries :        6 : Total  Size=        619 bytes  File Size  =        123 *
*Baskets :        1 : Basket Size=      32000 bytes  Compression=   1.05     *
*............................................................................*
*Br    2 :track_charge : track_charge/I                                      *
*Entries :        6 : Total  Size=        612 bytes  File Size  =        110 *
*Baskets :        1 : Basket Size=      32000 bytes  Compression=   1.00     *
*............................................................................*
root [3] myFirstTree->GetEntries()                                            → Treeに入っているentry数
(Long64_t) 6
root [4] myFirstTree->Show(0)                                                 → 0 entry目はなんですか?
======> EVENT:0
 track_pt        = 3200.43
 track_z         = 7643
 track_charge    = -1
root [5] myFirstTree->Scan()                                                 → entry順番にスキャンしてみます
************************************************
*    Row   *  track_pt *   track_z * track_cha *
************************************************
*        0 *   3200.43 *      7643 *        -1 *
*        1 *   2893.22 *      9834 *         1 *
*        2 *    3603.9 *     11232 *         1 *
*        3 *   9899.14 *     10232 *        -1 *
*        4 *   5674.32 *      8092 *         1 *
*        5 *   3432.33 *      7662 *         1 *
************************************************
(Long64_t) 6
root [6] myFirstTree->Draw("track_z")                                       → track_zという変数をヒストグラムにしてみます

簡単な確認なら上記でOKですが、実用的なことをする場合、最低でもマクロを書いてイベントループで解析する必要が出てきます。

#include <TFile.h>
#include <TTree.h>
#include <iostream>

void readNtuple()
{
  TFile* file = TFile::Open("myFirstTree.root", "READ");
  TTree* tree = dynamic_cast<TTree*>(file->Get("myFirstTree"));

  double track_pt;
  double track_z;
  int track_charge;

  tree->SetBranchAddress("track_pt", &track_pt);
  tree->SetBranchAddress("track_z", &track_z);
  tree->SetBranchAddress("track_charge", &track_charge);

  Long64_t nevents = tree->GetEntries();
  for(int i=0; i<nevents; i++) {
    tree->GetEvent(i);
    std::cout << i << ": pt=" << track_pt
          << ", z=" << track_z << ", charge=" << track_charge << std::endl;
  }
}

Exercise 4: Ntuple anaysis

2015年のATLASの1ランをL1TGCNtupleと言われる形式にして保存したものが
/home/junpei/ROOTtut2015/Exercise4/data15_13TeV.00280423.physics_Main.recon.NTUP.1.f629.00-00-18_L1TGCNtuple.Derivated.root
にあります。これを使ってミューオンのLevel-1 トリガー効率を出すことを目標とします。

L1TGCNtupleにどういう変数があるかなどはこちら. ただし今日用意したものはサイズの関係上、変数を一部落としています。

本当は Z粒子などを用いたTag-And-Probe法でバイアスが無くなるような解析が必要なんですが、 そこまでやる時間が今日はないのと、目的は実は違うところにあるので今回はパスします。
(なので結果は物理的には正しくはありません)

使う変数の説明、どうやるかはホワイトボードでも使って話すとして、 学んでほしいことは

  • TTree::MakeClass()
  • TTree::SetBranchStatus(...)
  • TH1::Sumw2()
  • TLorentzVector
  • efficiencyを計算するときにはどうするのか
  • TH1::Divide(...)
  • binominal error

(追記) 使う変数は

  • mu_n
  • mu_pt
  • mu_eta
  • mu_phi
  • trig_L1_mu_eta
  • trig_L1_mu_phi
  • trig_L1_mu_thrNumber

MakeClassで作ったマクロの使い方

     root [0] .L physics.C
     root [1] physics t
     root [2] t.Loop();

Reference


-- JumpeiMaeda - 2016-03-23

Edit | Attach | Watch | Print version | History: r16 < r15 < r14 < r13 < r12 | Backlinks | Raw View | WYSIWYG | More topic actions
Topic revision: r16 - 2016-03-30 - JumpeiMaeda
 
    • Cern Search Icon Cern Search
    • TWiki Search Icon TWiki Search
    • Google Search Icon Google Search

    Main All webs login

This site is powered by the TWiki collaboration platform Powered by PerlCopyright &© 2008-2019 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding TWiki? Send feedback