Page for understanding ANAROOT. Sorry for Japanese only.
Contents
はじめに
ここではANAROOTの初歩的な扱いに慣れた人がANAROOTの構造を理解するためのtutorial。 試験的に作成している段階なので変だと思う所は直接コメントを書いたりメールしてください。
演習1: ANAROOTをライブラリとして使う
ANAROOTは解析するためのライブラリ群で、これをコード上で呼ぶことで解析する。 ここではbookだとかstartだとかのANAROOTのコマンドは一切使わない。
step 1.0: 下準備
まず少しだけ前準備をしておく。 ただし、理研や東工大の解析マシンで
anarootlogon
などのshell commandがある場合は特に準備する必要はないはず。
ディレクトリ構造
ここでは解析を行うベースディレクトリを"mayoi"としたとき、
mayoi mayoi/db mayoi/macros mayoi/matrix mayoi/ridf mayoi/rootlogon.C mayoi/src
というディレクトリ構造を想定する。 公式に配布されているANAROOTの書庫をmayoi/src以下に展開し、make & make install すると ライブラリがmayoi/src/lib以下にinstallされる。 compile・install周りの詳しいことは公式サイトを見てください。
rootlogon.C
ANAROOTは単なるライブラリ群(libanaroot.so等, 動的ライブラリで検索)なので、 ROOTではこれを明示的に指定・呼び出ししてあげないと使いようがない。 そのために、例えばROOT起動時に実行される rootlogon.C にライブラリの場所の指定・呼び出しを書いといてあげる。
以下は愚直に書く場合の例。
void rootlogon()
{
gSystem->AddIncludePath("-Isrc/include");
gSystem->Load("libXMLParser.so");
gSystem->Load("src/lib/libanaroot.so");
}これをmayoi/rootlogon.Cとして保存しておく。 場合によってはカレントディレクトリにあるrootlogon.Cを読むように指定してない場合もあり、 その場合は.rootrcを編集する必要があるかもしれません。
step 1.1: 最小のプログラム
まずは最小のプログラムを理解しよう。 こちら(ArtAnalyzeDALI1.C)をダウンロードください。 このstep 1.1が本質で、あとは拡張するだけなので。
念のため実行方法を説明しておく。 cint上(rootを実行したときのコマンド入力画面)で実行する場合には以下の四つの方法がある。
// cint上でコンパイルしない方法 root[0] .x ArtAnalyzeDALI1.C // cint上でコンパイルする方法 root[0] .x ArtAnalyzeDALI1.C+ // ロードして実行 root[0] .L ArtAnalyzeDALI1.C root[1] ArtAnalyzeDALI1() // コンパイル&ロードして実行 root[0] .L ArtAnalyzeDALI1.C+ root[1] ArtAnalyzeDALI1()
cintは完全にc・c++の文法を理解できるわけではない(関数ポインタ、演算子のオーバーロード、テンプレート引数にオブジェクト等)のと自分の書いたコードのDEBUGという観点からコンパイルして実行することを勧める。 つまり、.x ArtAnalyzeDALI1.C+と最後に"+"をつけて実行する方法が推奨される。
では、説明を加えると、
#ifndef __CINT__ ... #endif
これは.x ArtAnalyzeDALI1.Cで呼んだときと、.x ArtAnalyzeDALI1.C+で呼んだときで ヘッダーファイルの読み込みがされない/されるを切り替えるため。 詳しくはプリプロセッサとかで調べてください。
"TArtHoge.hh"
TArtはANAROOT独自のオブジェクトということを表している単なる命名規則の一つ。 作法としてTArtHoge.hhにはTArtHogeクラスが定義されている (たまに継承したクラスも一緒にかかれていることがあるので注意)。 たとえばDALIのEnergyの情報を取ってきたいと思ったときに見るべきファイルはTArtCalibDALI.hhやTArtCalibDALI.ccとなる。 物理に関わるところは全てReconstructフォルダ以下にまとまっているが、見つからない時は、 ファイル検索であればfindやlocate、ファイルの中身の検索であればgrep等のlinux標準のコマンドを使えばよい。
TArtEventStore.hh: ridfファイルを読み込んでデコードしてrawdata型に詰めてくれる一番偉い人 TArtDALIParameters.hh: DALIのparameterを読み込んでDALIのparameter型に詰め込む TArtCalibDALI.hh: DALIへの窓口&本体。DALIという検出器その物を表したクラスではないということに注意。 TArtDALINaI.hh: DALIの素要素であるNaI一つのデータ形式 TArtDALIPara.hh: DALIのNaI用のparameterのデータ形式
ここで人によっては、DALI全体を表すクラスが独立にないことを気持ち悪いと思うかもしれないが、 諸事情により、素要素を定義するクラスの集合(vector<TArtDALINaI>の様なイメージ)を一つの検出器として扱っている。 具体的にはTClonesArrayというROOT標準のクラスに詰め込んでいる。 また一方ではbeam lineのplasticをTClonesArrayにまとめて突っ込んでたりする。 おそらくこの時点では何を言っているかわからないと思うが、 とりあえず構造がオブジェクト指向のセオリーに完全には乗っていないことだけは覚えていてほしい。 あと、Calibなんたらとついているものはあくまで”窓口兼本体”であることをよく頭に入れておくこと。
注意点をいくつか述べておく。
EventStore周りは深追いしても何も良い事はないので関数の中身を見て雰囲気を理解するぐらいにしておくべき。逆に、TArtCalibDALIだとかDALINaIとかはその先のコードまで理解できていることが望ましい。
- クラスをまたがるコンストラクタや関数内部で結構やっているので、その辺はとりあえず慣れてください。例えばTArtCalibDALIより先にTArtDALIParametersをコンストラクトしなければならない(コンストラクトする順番に依存するのは設計としては良くないとか突っ込まない)とか。
step 1.2: Histogramの定義
step 1.1にhistogramの定義を追加する。 こちら(ArtAnalyzeDALI2.C)をダウンロードください。
このマクロを回したら、.lsかls() すればhistogramが出来ているのが見えるはずである。 あとはhid->Draw()やhn()とすればちゃん表示されるはず。 (空っぽに見えるときはhenergyの定義域を4000にしたりlgyしたりridfファイル変えたり、 それでもだめならDEBUGと書いてあるところのコメントアウトを外してみてください)
多くのユーザーがちゃんとc++を勉強していない感があるので"new"について少し。 普通に変数を宣言してオブジェクトを生成する場合はスコープを抜けるときにそのオブジェクトは勝手に破棄(デストラクト)される。 一方newでオブジェクトを生成(コンストラクト)した場合はスコープを抜けてもそのオブジェクトは破棄されずメモリ上に残る。 例えばヒストグラム(TH1)はmacroで処理した後も参照したいのでメモリ上に残しておく必要があり、 敢えてnewして且つ敢えてdeleteしない。 これとは別にポインタの方が汎用的に扱えるという理由からnewしてしまう場合もあるが、 明確に目的が違うので混同しないこと。 TArtAnalyzeDALI2でnewしたTH1DのオブジェクトをTH1という親の型のポインタで受けているのは両方の目的がある。
もう一つ、ROOT標準の型の調べ方について。 基本的に"TH1 ROOT"とかでググれば公式に自動生成されているクラスリファレンスがhitする。 このページで、引数、関数名、戻り値、でそれっぽい関数はだいたい見つかる。 でも、例えばオブジェクトの名前を変える関数はTH1には載っていなくて、 これは親の関数の処理だから。 継承関係はClass ChartsのInheritanceに図で示されている。 子は親の処理はなんでも使えて、例えば名前の変更だったらTNamedというクラスにSetNameという名前で実装されている。 これを表すのにTNamed::SetNameって書くこともよくある。 継承とかはぱっと理解できるものではないので徐々に理解していってください。 動物である人間は動物ができることは人間もできるみたいなそういうごく自然な関係をプログラムで実現しているだけです。
あれやりたい、これやりたい、って漠然としたことは英語で"histogram title ROOT cern"とかで検索すればRoot Talkのスレッドが引っかかることも。
演習2: CalibHOGE周り
ここではCalibなんたらをいじっていきます。
上で使ったDALIのソースの中で計算すれば確かに何でも計算できるが、 これではみんながみんな似たような事を何回も書かなくてはいけなくなり、バグが組み合わせ的に増えていく(面倒という以上に問題)。 設計思想としては誰でも・いつでも必要とするものはTArtCalibDALIなどのクラスの中で計算しようということになっている。 誰でも・いつでもではないけどこの実験では、みたいなものはTAlなんたらの中で計算することがあるが、これは後ほど。 今だけちょこっと計算したいものは上記のようなmacroの中でちょろちょろっと計算する。 TArtCalib何たらはそのクラス単独で計算が簡潔するクラスで、 TArtReco何たらは他のCalibなどを使って計算する(粒子の運動量とか)クラス。 本質的な違いはないが実装の詳細は結構違うのでCalibに慣れてからRecoなんたらを考えるようにしてください。
step 2.0: 下準備
共有のsrcを使っている場合はそれを勝手に書き換えるわけにもいかないのでまずは下準備。 東工大の人向け。
共有のsrcを使っている人はおそらく
$ ls -la
すると、
src -> /home/samurai/exp/sm_com/anaroot/src
とかそんな感じでリンク先が表示されるかもしれません。 このリンクをやめて、実体に置き換えます。
$ rm src $ cp -r /svn/anaroot/branches/2012SAMURAI.dayone.tanaka ./src $ cd src $ ./autogen.sh --prefix=$PWD $ make $ make install
とかしてください。 ただし、かならず何をしているか考えた上でここに書いてあることを実行してください。
step 2.1: TArtCalibHOGEについて
ここでは、TArtCalibDALIを見ていきます。 まずはsource/Reconstruct/DALI/include/TArtCalibDALI.hhを開きましょう。 他のTArtCalibなんたらと比べるとわかりますが、 LoadData(), ClearData(), ReconstructData() は必ず実装されている関数です。 基本的にはrawdataをloadしてreconstruct(要は計算)してclearする、という単純な流れです。 つまり基本的にはReconstructData()を適当にいじればよいということ。 ではsource/Reconstruct/DALI/include/TArtCalibDALI.ccを開いてTArtCalibDALI::ReconstructDataを見てみましょう。
void TArtCalibDALI::ReconstructData() { // call after the raw data are loaded // このコメントはもう古いので無視
if(!fDataLoaded) LoadData(); // DALIのrawデータをもらってくる。RIDFからrawdataそのものを読み込むわけではない。
// 以下のまとまりについて解説すると、Tref(DALIの使っているTDCはジッターを除くためにTiming Referenceが必要)は
必ずfNaIArrayの最後に詰められる(ほんとにこの動作が保証されているかは知らない)ことを前提にして、
一番最後のNaIを取ってきてそこに格納されているTDC情報をtrefとして取ってきている。DALIに特有な処理なのでまぁあまり気にする必要はない。
// ----------------------------------------------
TArtDALINaI *nainai = 0; // DALIの要素であるNaI一個。ローカル変数の命名規則的にはTArtを除いて全て小文字にしたdalinaiとするべき。"= 0"はポインタの初期化。
Int_t tref = 0; // DALIの使っているTDCはジッターを除くためにTiming Referenceが必要(Decoderで差っ引くこともある)。
if(GetNumNaI() > 0){ // そのeventにおいて情報のあったNaIの数がGetNumNaIで返ってくる。情報が何もなかった場合のエラー処理。
nainai = (TArtDALINaI*)fNaIArray->At(GetNumNaI()-1); // fNaIArrayはそのeventにおいて情報のあったNaIが詰め込まれた配列みたいなもの。TClonesArrayというROOT標準のクラス。かなり特殊なのでROOTに慣れていない人はとりあえずvectorのようなものだと思っておけばよい。At()はその順番(IDではなくfNaIArrayに詰められた順番)のNaIを返す。戻り値の型はTObject型なのでTArtDALINaI型に明示的にキャストしてる。これはダウンキャストと呼ばれる危険なキャストなので知識がなければ放置。
tref = nainai->GetRawTDC(); // さっき取ってきたNaIに格納されているTDCのrawデータを取り出している。
}
// ----------------------------------------------
for(Int_t i=0;i<GetNumNaI();i++){ // 情報のあった全てのNaIに対してloop処理
TArtDALINaI *nai = (TArtDALINaI*)fNaIArray->At(i); // i番目に詰め込まれたNaIを取ってくる。
TArtDALINaIPara *para = fNaIParaArray[i]; // i番目に詰め込まれたNaIに対応するParameterを取ってくる。
Int_t adc = nai->GetRawADC(); // ADCのデータを取ってくる
Int_t tdc = nai->GetRawTDC(); // TDCのデータを取ってくる
Double_t theta = nai->GetTheta(); // parameterから取ってこれるがNaIにコピーされたthetaを使っている。個人的にはどうかと思う。
Double_t costheta = TMath::Cos(theta*3.1415926535/180.); // cos(theta)をよく使うので。
Double_t fEnergy = adc * para->GetQCal() + para->GetQPed(); // Energyキャリブレーション
nai->SetEnergyWithoutT(fEnergy); // T情報があってもなくてもenergy情報を格納する。pedestalとか？
nai->SetCosTheta(costheta); // 計算したcos(theta)を格納
if(tdc>0 && adc>0 && tdc<100000 && adc<4000){ } // TDC情報とADC情報がまっとうなら特に何もしない。書き方意味不明。ADCが4000で切っているのも意味不明。
else{ continue; } // 上の条件以外だと次のloopへ。
nai->SetEnergy(fEnergy); // energyを格納
Double_t fTime = (double)tdc * para->GetTCal(); // TCal
Double_t fTimeOffseted = fTime + para->GetTOffset(); // 時間のoffset調整
// 後はデータを格納して終わり
nai->SetTime(fTime);
nai->SetTimeOffseted(fTimeOffseted);
if(fTimeOffseted > 180. && fTimeOffseted < 240.){
nai->SetTimeTrueEnergy(fEnergy);
nai->SetTimeTrueTime(fTime);
nai->SetTimeTrueTimeOffseted(fTimeOffseted);
}else{
nai->SetTimeTrueEnergy(-1.);
nai->SetTimeTrueTime(-1000.);
nai->SetTimeTrueTimeOffseted(-1000.);
}
}
fReconstructed = true; // ANAROOTのお約束としてフラグ立てしておく。
return;
}ということでした。
step 2.2: 新しい物理量の追加
TArtCalibDALIの中で既に計算・格納されている物理量(energyとか)の計算方法を変えるとかはいいとして、 ここでは新しい物理量を計算して格納することをやってみる。 具体的にはbetaをとりあえず0.6として、Reconstructのなかでドップラー補正してみる。
まずTArtDALINaI.hhに新しく追加したい物理量のためのメンバ変数を追加する。 追加するのはドップラー補正後のenergyだが、既にfDoppCorEnergyという変数があるので、 とりあえずfDoppCorEnergy2という名前にしよう。 TArtDALINaI.hhの変更点だが、簡単にはfDoppCorEnergyを検索してfDoppCorEnergy用に定義されている関数や 初期化処理等を真似ればよい。
TArtDALINaI() : ..., fDoppCorEnergy2(-1), ...
virtual void SetDoppCorEnergy2(Double_t val){fDoppCorEnergy2 = val;}
virtual Double_t GetDoppCorEnergy2(){return fDoppCorEnergy2;}
Double_t fDoppCorEnergy2;つぎにTArtCalibDALI.ccを開いて
void TArtCalibDALI::ReconstructData()
{
...
for(Int_t i=0;i<GetNumNaI();i++){
...
Double_t beta = 0.6;
Double_t doppcorenergy2 = fEnergy*(1-beta*costheta)/sqrt(1-beta*beta);
nai->SetDoppCorEnergy2(doppcorenergy2);
...
}
}と適当にかく(ドップラー補正がこれでいいのかは知らない)。 あとはsrcディレクトリに戻ってmake, make installする。 これでmacroにnai->GetDoppCorEnergy2()とか書けばReconstructの中で計算したものが見えるはず。 ただし、ヘッダファイルの変更はmakeしてもその更新が反映されないことがある (たぶんROOT用のディクショナリに反映されていないだけな気がする)。 こういうときはmake cleanしてみるといいかも。 それでもダメなときは違うsrcを編集していたってこともあったり。 マクロをコンパイルして使っている場合はそのコンパイルしてできるhoge_C.dとhoge_C.soを消すといい場合も。 実行時エラーで落ちた場合にhoge_C.dとhoge_C.soが変に残って変更が反映されないとかたまにあります。
...ということでできたでしょうか。
step 2.3: parameterの扱い
上でbetaを直接0.6と書いたが、これだとあんまりなので、parameter fileから呼んでみる。
step 2.4: 新しいクラスの追加
新しいクラスの追加。 適当に真似すればなんとかなる。 ただし、anahoge_linkdef.hhにちょこっと書く必要があるのを忘れずに(ROOTの辞書への登録)。
番外編: svn
上記でソースを変更するためにコピったと思うが、これだと本流のソースが更新した時にその更新を反映しにくい。 一応大体のバグはとれているので本流の更新を反映しなくてもやっていけるが、 例えばDCのtrackingの方法が改善したとか、そういうのをきっと自前のソースに反映したくなる。 ANAROOTはsvnという仕組みで管理されていて、svnのコマンドを使うと 手元のソースをbranch(支流)としてサーバーにuploadしたり、本流をbranchに反映したりできる。 ただし使い方の思想にかなり癖があり、初めてだとさっぱりということもしばしば (不便というわけではない。むしろものすごく便利)。 必要なら自分で頑張って理解するか誰かに聞いてください。
演習3: TAl(読み方: ティーアナループ, タル)の使い方
TAlのいじり方。TAlはTArtAnaLoopの略称。 TAlDALIExampleとTAlEncDALIExampleがあるが、 Enc(フランス人は”エンク”ではなく"エ"と"オ"の中間音で発音する)とついたものは 特別にanafileで扱えるようにしたクラス。 そうでないTAlDALIExampleとかはanafileは使えないがbook,start,stop等の関数で扱える。
step 3.1: MacroでTAlDALIExample
bookで読み込めばANAROOTに実装されているstartやstop,status等の関数が使えるので ちょこっと便利になる。
step 3.2: MacroでTAlEncDALIExample
これはほとんど趣味の世界。 srcを変更せずにanafileで扱える物理量やanalysを増やせる。
演習4: Treeを作る・解析する
常にRaw Dataから出発して解析しているとものすごく解析に時間が掛かる。 そこでROOT標準のTreeを用いる。 Treeはかなり高性能だが、ここでの使い道はANAPAWユーザー的にはDumpファイル?に相当する。 Treeをどこでどう作るかというのはいくつかやりようがあるので紹介する。
step 4.1: Macroで直に
まぁ好きにやってください。
step 4.2: TAlに仕込む
analoopの項に書いた気がする。 適当にやるならConstructの中でnewするだけ。 真面目にやるならbookをオーバーロードしてその中でnewする。 2回以上bookする場合に違いが出てくるだけなのでofflineではConstructのなかで treeを作ってしまえばよいかも。 ただし、TAlの中でROOT標準のWrite関数は読んではだめ。フリーズします。
step 4.3: anafileを使って
2012/08/05現在、2012SAMURAI.dayone.tanakaでのみanafileでtreeが作れるようになってます。 anafileに、
// ana/nebula.ana <analys> 15, <branch> 101, 15,1,144, 1, "id" "nebulatree" <branch> 101, 15,1,144,10, "quraw" "nebulatree" <branch> 101, 15,0, 0,51, "multiplicity" "nebulatree" <filltree> "nebulatree"
みたいな感じに書くだけで簡単にtreeが作れます。 idが範囲指定されているときは自動でvector<Dobule_t>になります。 範囲が1のときはDouble_tになります。 結構便利です。 Calibrationでアイテレーションしたいとかslew補正をいろいろ試したいとか、そういうときには重宝します。 また、物理量が出るまで行ったときでも並行してヒストグラムを定義してcheckできるので便利かと思います。
特定のクラスをボッコリ詰めるのにはまだ対応してませんが暇があれば実装するかもしれません(できるのかな...)。
step 4.4: treeの解析
まぁ。