PAW
メモ
- bbtcal2.kumac のフィットは二次でフィットしているので、気を付ける。
- データを線形フィットする
v/cr x(3) r 1.0 2.0 3.0 v/cr y(3) r 2.1 4.1 5.9 v/f x y ! p1
- データ点を * でプロットするには以下のように、* オプションをつける
v/f x y ! p1 *
- データを線形フィットし内挿(外挿)する
v/cr x(3) r 1.0 2.0 3.0 v/cr y(3) r 2.1 4.1 5.9 v/cr p(2) v/f x y ! p1 * 2 p * x=1.5 のときの y の値を内挿 mess $sigma(p(1)+p(2)*1.5) * x=4.0 のときの y の値を外挿 mess $sigma(p(1)+p(2)*4.0) * y=3.5 のときの x の値を逆関数から内挿 mess $sigma((3.5-p(1))/p(2))
- データを任意関数でフィットする
macro fit exe def_func v/cr x(3) r 1.0 2.0 3.0 v/cr y(3) r 2.1 2.9 3.8 v/cr p(2) r 0.0 1.0 v/f x y ! fit_func * 2 p return macro def_func app comis quit real function fit_func(x) common/pawpar/par(2) real x fit_func=par(1)+par(2)*x end quit return
- sigma を用いた vector の計算では、8文字以上の名前の vector は扱えない
- sigma x=x/vsum(x) とかもOK
- kumac で良く使うコード
xbins = $hinfo([hid],'xbins') xmax = $hinfo([hid],'xmax') xmin = $hinfo([hid],'xmin') dx = ([xmax] - [xmin])/[xbins]
- kumac で使うコード
macro xwinsize for str in xsiz xmgl xmgr xwin xzones [str]=$grafinfo([str]) endfor xwinsize = [xsiz] - ([xmgl] + [xmgr] + [xwin]*[xzones]) return [xwinsize] macro ywinsize for str in ysiz ymgu ymgl ywin yzones [str]=$grafinfo([str]) endfor ywinsize = [ysiz] - ([ymgu] + [ymgl] + [ywin]*[yzones]) return [ywinsize]
- COMIS 内では、コメント ! は必ず、文の後に書く。
- COMIS で定義した関数に配列を渡す
vec/cre x(2) r 3 4 app comis quit real function arr_sum(arr,narr) integer narr real arr(narr) arr_sum=arr(2)+arr(1) end program main vector x write(*,*) arr_sum(x,2) end quit
- kumac 内での文字列。
str1 = aa str2 = str1 str3 = 'bb' str4 = 'str3' str5 = [str1] str6 = [[str2]] str7 = '[str1]' str8 = '[[str2]]' str9 = [str3] str10 = [[str4]] str11 = '[str3]' str12 = '[[str4]]' str13 = $quote([str5]) str14 = $quote([str6]) str15 = $quote([str7]) str16 = $quote([str8]) str17 = $quote([str9]) str18 = $quote([str10]) str19 = $quote([str11]) str20 = $quote([str12]) str21 = $quote([str1]) str22 = $quote([[str2]]) str23 = $quote('[str1]') str24 = $quote('[[str2]]') str25 = $quote([str3]) str26 = $quote([[str4]]) str27 = $quote('[str3]') str28 = $quote('[[str4]]') mess aa [str1] $quote([str1]) | aa aa aa と表示される mess str1 [str2] $quote([str2]) | str1 str1 str1 と表示される mess 'bb' [str3] $quote([str3]) | bb bb bb と表示される mess 'str3' [str4] $quote([str4]) | str3 str3 str3 と表示される mess [str1] [str5] $quote([str5]) | aa aa aa と表示される mess [[str2]] [str6] $quote([str6]) | aa aa aa と表示される mess '[str1]' [str7] $quote([str7]) | [str1] aa [str1] と表示される mess '[[str2]]' [str8] $quote([str8]) | [[str2]] aa [[str2]] と表示される mess [str3] [str9] $quote([str9]) | bb bb bb と表示される mess [[str4]] [str10] $quote([str10]) | bb bb bb と表示される mess '[str3]' [str11] $quote([str11]) | [str3] bb [str3] と表示される mess '[[str4]]' [str12] $quote([str12]) | [[str4]] bb [[str4]] と表示される mess $quote([str5]) [str13] $quote([str13]) | aa aa aa と表示される mess $quote([str6]) [str14] $quote([str14]) | aa aa aa と表示される mess $quote([str7]) [str15] $quote([str15]) | [str1] [str1] [str1] と表示される mess $quote([str8]) [str16] $quote([str16]) | [[str2]] [[str2]] [[str2]] と表示される mess $quote([str9]) [str17] $quote([str17]) | bb bb bb と表示される mess $quote([str10]) [str18] $quote([str18]) | bb bb bb と表示される mess $quote([str11]) [str19] $quote([str19]) | [str3] [str3] [str3] と表示される mess $quote([str12]) [str20] $quote([str20]) | [[str4]] [[str4]] [[str4]] と表示される mess $quote([str1]) [str21] $quote([str21]) | aa aa aa と表示される mess $quote([[str2]]) [str22] $quote([str22]) | aa aa aa と表示される mess $quote('[str1]') [str23] $quote([str23]) | [str1] [str1] [str1] と表示される mess $quote('[[str2]]') [str24] $quote([str24]) | [[str2]] [[str2]] [[str2]] と表示される mess $quote([str3]) [str25] $quote([str25]) | bb bb bb と表示される mess $quote([[str4]]) [str26] $quote([str26]) | bb bb bb と表示される mess $quote('[str3]') [str27] $quote([str27]) | [str3] [str3] [str3] と表示される mess $quote('[[str4]]') [str28] $quote([str28]) | [[str4]] [[str4]] [[str4]] と表示される
- str1 = aa, str3 = 'bb' は同じ意味になるようだ。代入時にクオートが外される?
- mess [str7] と mess [str23] で, 結果が aa と '[str1]' と違うのはなんで?代入時にクオートが外されるから?
- ネストされた参照の出力。
str1 = aa str2 = str1 str3 = bb str4 = [[str2]] str5 = [[str3]] mess [[str2]] [[str3]] [str4] [str5] $quote([str4]) $quote([str5]) * aa [bb] aa [bb] aa [[str3]] と表示される
- 未定義の文字列がネスト内にあると、参照自体が行われない。ただし、 mess 内は特別?
- 未定義の変数の判別。
str1 = aa str2 = str1 str3 = bb str4 = [[str2]] str5 = [[str3]] mess [str4] '[[str2]]' if ([str4] .eq. '[[str2]]') then mess '[str4] eq!' else mess '[str4] not eq!' endif mess [str5] '[[str3]]' if ([str5] .eq. '[[str3]]') then mess '[str5] eq!' else mess '[str5] not eq!' endif 実行結果 aa [[str2]] [str4] not eq! [bb] [[str3]] [str5] eq!
- データを書き出し、COMIS で読み込むマクロ
macro main APPlication DATA dat.txt 2.2 3.3 4.4 dat.txt APPLication COMIS QUIT program demo integer i,j real aa(10) i=1 open(unit=60,file='dat.txt',status='old') do while 1 read(60,*,end=30) aa(i) i=i+1 enddo 30 continue close(60) do j=1,i-1 write(*,*) j,aa(j) enddo end QUIT return
- vec/fit の引数
- 以下の感じでもいける。ey を ! としている。
vec/cre x(3) r 1 2 3 vec/cre y(3,2) r 2 3 4 5 6 7 vec/fit x y(1:3,2) ! p1
- COMIS で call b('abcdef'//'ghijkl') という呼びだしはできない。
- Minuit の chisquare
- CERN Program Library Long Writeup D506
- MINUIT - Function Minimization and Error Analysis Reference Manual Version 94.1 p.4 (1.1)式
: カイ二乗
: フィットのパラメータ
: フィットする関数
: フィットしたいデータ
: フィットしたいデータのエラー
- PAW 画面の再描画
window が重なってもヒストグラムがまた表示されるようにする。参考 : PAW misc Tips
/etc/X11/xorg.conf の Section "Screen" 〜 EndSection の間に、以下の行を追加し、X を再起動する。(ログアウトしてログインすれば、X が再起動されるかも?)
Option "backingstore"
- SL5.5 (32 bit) だと、ワークスペースを変えて戻ってもちゃんと表示されるが、SL6.4 (64 bit) だとワークスペースを変えて戻ると画像が消える。なんで?
kumac 中の整数は 231-1 = 2147483647 まで表示可能
- X11-interface-only の IGSET value
(参考 : CERNLIB documentation - HIGZ/HPLOT - The graphic macroprimitives - Setting attributes)
- DRMD : Drawing mode: 1.=copy 2.=xor 3.=invert
- SYNC : Synchronise the graphics in X11 1.=yes 0.=no
- 2BUF : 10*(WKID)+(double buffer mode: 1.=on 0.=off)
- WKID : Workstation identifier
データをプロット & 3次関数でフィット & 途中の値を表示 するマクロ
macro sub sigma x=array(10,250.#2500.) vec/cre y(10) r 50.2 32.4 24.0 19.2 16.1 14.1 12.5 11.2 10.0 9.33 * graph 10 x y vec/cre par(4) r vec/fit x(3:7) y(3:7) ! p3 ! 0 par fun1 101 55.251-0.63695E-01*x+0.34571E-04*x**2-0.69333E-08*x**3 100 0. 3000. s sigma x=1173 mess $sigma(55.251-0.63695E-01*x+0.34571E-04*x**2-0.69333E-08*x**3) sigma x=1332 mess $sigma(55.251-0.63695E-01*x+0.34571E-04*x**2-0.69333E-08*x**3) return
- FIT, SET, OPT の端末への出力をファイルに書き出す
PAW> output 60 output.dat ! output.dat というファイルに set コマンドと opt コマンドの結果を書き出す PAW> set PAW> opt PAW> output -1 ! ファイルを閉じて、出力先を元に戻す
- kumac でやりたいこと
- 以下の様なことがやりたいけどできない。[v] を $unquote([v]) にしてもダメ。
v='aa bb cc' for str in [v] mess [str] endfor
- 以下のことはできる。[*] は [1] [2] [3] [4] ...
for str in [*] mess [str] endfor
疑問
- histogram に graph コマンドで上書きできるか?
- chbin て結局なにやってる?(ANAPAW)
- ライブラリをリンクできる?
データ変換
似たような図として、CERN Program Library Long Writeup Q121 PAW Physics Analysis Workstation An Introductory Tutorial p.28の Figure 2.2: PAW entities and their related commands がある。
① vector -> histogram
vector データを histogram に変換するには、 VECTOR/DRAW コマンドを用いる。
* sample_vect という名前の real 型 vector (要素数は5) を作り、初期値として 1 2 3 2 1 を代入 vec/cre sample_vect(5) r 1 2 3 2 1 * vector のデータをヒストグラムに変換。x軸の bin 数、最大値、最小値は自動的に決まる。 vec/draw sample_vect 100 * 試しに histogram をプロット。 hi/pl 100
ただし、VECTOR/DRAW コマンドでは、x軸の bin 数、最大値、最小値は自由に決められないかも?自分で決めたいときは、HISTOGRAM/PUT_VECT/CONTENTS コマンドを用いる。
* sample_vect という名前の real 型 vector (要素数は5) を作り、初期値として 1 2 3 2 1 を代入 vec/cre sample_vect(5) r 1 2 3 2 1 * sample_hist という名前の 1次元 histogram を作る。id を101、x軸の bin 数を5、最小値を0、最大値を10とする。 1dhisto 101 'sample_hist' 5 0 10 * 上で作った histogram (id=101) を表示するとき、y軸の表示範囲を指定したい場合、min と max コマンドで指定する。以下では、y軸の最小値を0.5、最大値を3.5 としている。指定しなくてもよい。 min 101 0.5 max 101 3.5 * vector のデータを histogram に変換。 put/con 101 sample_vect * 試しに histogram をプロット。 hi/pl 101
② histogram -> vector
histogram の各 bin の値を vector に代入するには、 HISTOGRAM/GET_VECT/CONTENTS コマンドを用いる。ここでは、 vector -> histgram のところで作った id=101 の histogram がすでにあるとする。
* histogram があるか確認。 hi/list * 以下のように表示されれば良い。 * ===> Directory : //PAWC * 101 (1) sample_hist * * sample_hist の bin 数を確認。 mess $hinfo(101,'xbins') * 5 と表示される。 * hist2vect という名前の real 型 vector (要素数は5) を作る。 vec/cre hist2vect(5) r * 上記コマンドの代わりに、vec/cre hist2vect($hinfo(101,'xbins')) r としてもよい。 * vector のデータをヒストグラムに変換。 get/con 101 hist2vect * vector の中身を確認。 vec/print hist2vect * 以下のように表示される。 * HIST2VECT(1) = 1 * HIST2VECT(2) = 2 * HIST2VECT(3) = 3 * HIST2VECT(4) = 2 * HIST2VECT(5) = 1
③ vector -> text
vector の中身を text ファイル (ASCII)に書き出すには、 VECTOR/WRITE コマンドを用いる。
* x, y という名前の real 型 vector (要素数は5) を作る。 vec/cre x(5) r 1 2 3 4 5 vec/cre y(5) r 1 4 9 16 25 * vector x, y の中身を xy.dat という名前の text ファイルに書き出す。書き出すフォーマットは、4桁中に x の データを小数点1桁まで書き出し、6桁中に y の データを小数点2桁まで書き出すようにする。さらに、xとyのデータの間にスペースを2桁挿入する(2Xの部分に相当)。フォーマットは省略不可か?省略可能かも? vec/write x,y xy.dat 'F4.1,2X,F6.2' * 書き出した text を見てみる。 sh cat xy.dat * 以下のように表示される。 * 1.0 1.00 * 2.0 4.00 * 3.0 9.00 * 4.0 16.00 * 5.0 25.00
④ text -> vector
text ファイル (ASCII)を読み込んで vector に代入するには、 VECTOR/READ コマンドを用いる。
* vector -> text で作った xy.dat というデータファイルがワーキングディレクトリにあるとする。 sh cat xy.dat * 以下のように表示される。 * 1.0 1.00 * 2.0 4.00 * 3.0 9.00 * 4.0 16.00 * 5.0 25.00 * データを入れるための、 x, y という名前の real 型 vector (要素数は5) を作る。 vec/cre x(5) r vec/cre y(5) r * xy.dat から、vector x, y に値を代入。'F4.1,2X,F6.2' の部分は読み込みのフォーマットを表すが、この場合は書かなくても良い。 vec/read x,y xy.dat 'F4.1,2X,F6.2' * vector x, y の中身を確認。 vec/print x * 以下のように表示される。 * X(1) = 1 * X(2) = 2 * X(3) = 3 * X(4) = 4 * X(5) = 5 vec/print y * Y(1) = 1 * Y(2) = 4 * Y(3) = 9 * Y(4) = 16 * Y(5) = 25
⑤ text -> ntuple
text ファイル (ASCII)を読み込んで ntuple に代入するには、 NTUPLE/READ コマンドを用いる。
* vector -> text で作った xy.dat というデータファイルがワーキングディレクトリにあるとする。 sh cat xy.dat * 以下のように表示される。 * 1.0 1.00 * 2.0 4.00 * 3.0 9.00 * 4.0 16.00 * 5.0 25.00 * データを入れるための ntuple を作る。 ntuple の id は102, 名前は sample_ntpl, データの種類は 2 (xとyで2つ), データ数は 5 である。 nt/cre 102 sample_ntpl 2 ' ' 5 xx yy * xy.dat から、ntuple に値を代入。 nt/read 102 xy.dat 'F4.1,2X,F6.2' * 'F4.1,2X,F6.2' の部分は読み込みのフォーマットを表すが、この場合は書かなくても良い。 * 試しに ntuple をプロットする。xx vs. yy の 2次元ヒストグラムをプロット。 nt/plot 102.yy%xx * yy の重みがついた xx の 1次元ヒストグラムをプロット。 nt/plot 102.xx yy
⑥ ntuple -> histogram
ntuple から histogram を作るには、 NTUPLE/PROJECT コマンドを用いる。
2次元 ntuple -> 2次元 histogram
* text -> ntuple で作った id=102 の ntuple があるとする。 nt/list * * ===> Directory : //PAWC * 102 (N) sample_ntpl * データを入れるための2次元 histogram を作る。 * histogram の id は103, 名前は sample_2dhisto, x, y軸の bin 数を30、最小値を0、最大値を30とする。 2dhist 103 sample_2dhist 30 0 30 30 0 30 * ntuple から histogram に値を代入。 nt/proj 103 102.yy%xx * 試しに histogram をプロットする。 hi/plot 103 colz
2次元 ntuple -> 1次元 histogram (重みありの histogram を作る)
* text -> ntuple で作った id=102 の ntuple があるとする。 nt/list * * ===> Directory : //PAWC * 102 (N) sample_ntpl * データを入れるための1次元 histogram を作る。 * histogram の id は104, 名前は sample_1dhisto, xの bin 数を5、最小値を0.5、最大値を5.5とする。 1dhist 104 sample_1dhist 5 0.5 5.5 * ntuple から histogram に値を代入。 nt/proj 104 102.xx yy * 試しに histogram をプロットする。 hi/plot 104
⑦ vector -> graph
vector から graph をプロットする場合、エラーを付けないときは GRAPHICS/PRIMITIVES/GRAPH を用い、エラーを付けるときは GRAPHICS/HPLOT/ERRORS を用いる。
- エラーバー無し (こっちの方が楽)
* x, y という名前の real 型 vector (要素数は5) を作る。 vec/cre x(5) r 1 2 3 4 5 vec/cre y(5) r 1 4 9 16 25 * プロット graph 5 x y * graph $vlen(x,1) x y でも良い。
- エラーバー有り (少し面倒)
* x, y という名前の real 型 vector (要素数は5) を作る。 * さらに、エラーの値を入れる vector として、ex, ey という名前の real 型 vector (要素数は5) を作る。 vec/cre x(5) r 1 2 3 4 5 vec/cre y(5) r 1 4 9 16 25 vec/cre ex(5) r 5*0.3 vec/cre ey(5) r 5*1 * プロット hp/err x y ex ey 5 ! ! W * hp/err x y ex ey $VLEN(x,1) ! ! W でも良い。
⑧ raw data (vector) -> histogram
vector 形式の raw データを histogram に詰めていくには、 VECTOR/HFILL コマンドを用いる。
* sample_vector という名前の real 型 vector (要素数は5) を作り、初期値として 1 2 3 2 1 を代入 vec/cre sample_vector(5) r 1 2 3 2 1 * vector のデータを入れるために、sample_histo という名前の histogram を用意する。histogram の ID は100, x軸のbin の数は3, 最大、最小値はそれぞれ 0.5, 3.5。 1dhisto 100 'sample_histo' 3 0.5 3.5 * vector のデータをヒストグラムに詰める。 vec/hfill sample_vector 100 * 試しに histogram をプロットする。 hi/pl 100