2014年11月11日火曜日

TSplineの境界条件の設定の方法

ROOTで点を補間しようと思った時に、点のある区間はスプライン補間で行い、その外は別の関数で補間することを考えた。
その際に、端点でスプラインと任意の関数が滑らかにつながるようにしたい。
このとき、TSplineに境界条件を入れることを考える。
CernのTSpline3の説明を見ると、

TSpline3(const char* title, const TGraph* g, const char* opt = 0, Double_t valbeg = 0, Double_t valend = 0)

となっており、valbegとvalendに端点の1階微分を境界条件として入れてあげれば良さそう。
デフォルトでは0となっているoptを設定しなければならない。
このoptに何を設定するかがわかりづらかった。
どうやらこのページを見ると、optにcharで"b1e1"と入れると、1階微分の境界条件が設定されるようだ。

試しにやってみたのが以下のプログラム(pythonで動かしてみた)。
ここでは、外の点は端から2点をつないだ直線でつなぐことを考えて、境界条件として端から2点の傾きを求めた。

>>> import ROOT as pR
>>> import array
>>> a = array.array("d",[1,2,3,4,5]) 
>>> b = array.array("d",[10,3,7,4,8])
>>> g = pR.TGraph(len(a),a,b)
>>> g.SetMarkerStyle(2)
>>> s = pR.TSpline3("s",g) # 境界条件を設定しないスプライン
>>> s.SetLineColor(2) #赤色
>>> tmp = [pR.Double(),pR.Double(),pR.Double(),pR.Double()]
>>> # 左端の傾き
>>> g.GetPoint(0,tmp[0],tmp[1]) # 始点
>>> g.GetPoint(1,tmp[2],tmp[3]) # 始点のお隣
>>> derF = (tmp[3]-tmp[1])/(tmp[2]-tmp[0])
>>> # 右端の傾き
>>> g.GetPoint(g.GetN()-2,tmp[0],tmp[1]) # 終点一つ前
>>> g.GetPoint(g.GetN()-1,tmp[2],tmp[3]) # 終点
>>> derE = (tmp[3]-tmp[1])/(tmp[2]-tmp[0])
>>> s2 = pR.TSpline3("s2",g,"b1e1",derF,derE) # 境界条件を設定したスプライン
>>> s2.SetLineColor(3) # 緑色
>>> 
>>> g.Draw("AP")
>>> s.Draw("same")
>>> s2.Draw("same")


2014年10月22日水曜日

gnuplot小技

gnuplotでの小技のまとめ

1, エラーバーの端っこ消す

set bar 0

2, 2つの関数で囲んだ部分をfillする

plot "+" us 1:(f(\$1)):(g(\$1)) with filled curves closed

注意すべきことは、関数f(x), g(x)を括弧でくくること。
そうしないと、undefined value といわれる。

2014年5月21日水曜日

PyROOTを使いこなすためには…?

PyROOTを使ってROOTファイルの解析を行う。

設定方法はGoogleで検索すれば見つかるでしょう…
問題は使い方!そんなに普及していない(?)のか日本語記事がなかなか見つからない!!
ROOT自体を使いこなしている訳ではないので,Python内で何のオブジェクトを使えばいいのかに結構悩んだ。。。

ってことでROOTを使い始めて,というかPyROOTを使い始めて数ヶ月経ったので,新しくROOTファイルを作る方法について重点をおいてまとめてみる。


・PythonにROOTをimport
>>> import PyROOT as pR

・ROOTファイルを開く
>>> # 既存のfileを開く
>>> f_in = pR.TFile("test.root")
>>> # 新しいroot fileを作る
>>> f_new = pR.TFile("output.root","recreate")

・1次元ヒストグラム
TH1を使う。
>>> # 0から1024の間を1024ビンに区切ったヒストグラム
>>> h1 = pR.TH1D("h1","h1",1024,0,1024)
>>> h1.SetBinContent(10,10)  # 10ビン目に10つめる 
>>> h1.Fill(350) # 350ビン目に1つだけつめる。
>>> h1.Draw()    # 描画
以上の操作で以下の図ができる。

・ディレクトリ構造を作る
>>> # カレントディレクトリの表示
>>> f_new.pwd()
output.root:/
>>>
>>>
>>> # ディレクトリ操作
>>> cur_dir = pR.gDirectory
>>> cur_dir.ls()
TFile**  output.root  
 TFile*  output.root   
>>> cur_dir.pwd()
output.root:/
>>>
>>>
>>> # ディレクトリ作成
>>> subdir = f_new.mkdir("subdir")
>>> f_new.ls()
TFile**  output.root  
 TFile*  output.root   
  TDirectoryFile*  subdir subdir  
  KEY: TDirectoryFile subdir;1 subdir
>>>
>>>
>>> # ディレクトリ移動
>>> f_new.cd("subdir")
True
>>> cur_dir.pwd()
output.root:/subdir
>>> f_new.pwd()
output.root:/
>>>
>>>
>>> # subdir 内に平均が512で分散が10となる
>>> # ガウス分布のヒストグラムを作る
>>> h1 = pR.TH1D("h1","h1",1024,0,1024)
>>> for i in range(1000) :
...    h1.Fill(pR.gRandom.Gaus(512,10)
>>> cur_dir.ls()
TDirectoryFile*       subdir subdir
 OBJ: TH1D      h1    h1 : 0 at: 0x7faa9aeec930 
>>>
>>>
>>> # 上位ディレクトリへ移動
>>> f_new.cd("../")
>>> cur_dir.pwd()
output.root:/
>>> f_new.ls()
TFile**         output.root
 TFile*         output.root
  TDirectoryFile*             subdir subdir
   OBJ: TH1D    h1    h1 : 0 at: 0x7faa9aeec930
>>>
>>>
>>> # 平均が512で分散が10となる
>>> # ガウス分布のヒストグラムを作る
>>> h2 = pR.TH1D("h2","h2",1024,0,1024)
>>> for i in range(1000) :
...    h2.Fill(pR.gRandom.Gaus(512,10)
>>> f_new.ls()
TFile**         output.root
 TFile*         output.root
  TDirectoryFile*             subdir subdir
   OBJ: TH1D    h1    h1 : 0 at: 0x7faa9aeec930
  OBJ: TH1D     h2    h2 : 0 at: 0x7faa9e10e730
>>>
>>>
>>> # ヒストグラムを作る
>>> h3 = pR.TH1D("h3","h3",1024,0,1024)
>>> f_new.ls()
TFile**         output.root
 TFile*         output.root
  TDirectoryFile*             subdir subdir
   OBJ: TH1D    h1    h1 : 0 at: 0x7faa9aeec930
  OBJ: TH1D     h2    h2 : 0 at: 0x7faa9e10e730 
  OBJ: TH1D     h3    h3 : 0 at: 0x7faa9dba1ec0
>>>
>>>
>>> # ヒストグラムの削除
>>> h3.Delete()
>>> f_new.ls()
TFile**         output.root
 TFile*         output.root
  TDirectoryFile*             subdir subdir
   OBJ: TH1D    h1    h1 : 0 at: 0x7faa9aeec930
  OBJ: TH1D     h2    h2 : 0 at: 0x7faa9e10e730
>>>
>>>
>>> # ヒストグラムをファイルへ書き込み,閉じる
>>> f_new.Write()
>>> f_new.Close()

てな感じかな。やり方にはいろいろあるだろうけど自分が試してみてできたやり方が上の方法。(なんかみにくいね…)