Locked History Actions

PAW

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 *

    fit.png

  • データを線形フィットし内挿(外挿)する
    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

    fit_func.png

  • 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)式
  • $$\chi^2(\alpha)=\sum_{i=1}^{n}\frac{(f(x_i,\alpha)-e_i)^2}{\sigma_i^2}$$

  • $$\chi^2$$ : カイ二乗

  • $$\alpha$$ : フィットのパラメータ

  • $$f(x_i,\alpha)$$ : フィットする関数

  • $$e_i$$ : フィットしたいデータ

  • $$\sigma_i$$ : フィットしたいデータのエラー

  • 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 の端末への出力をファイルに書き出す
    • 参考 : How to print the output of HIST/FIT into a text file ?

    • 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)
  • ライブラリをリンクできる?

データ変換

paw.png

似たような図として、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