モンテカルロ積分(N次元球体積)

Steering fileを自分で書いてみて実際にArtemisを動かしてみよう。

チュートリアルではこの記事で準備した解析環境を使用します。

ここを参考に、4次元球までの体積をモンテカルロ積分で計算してみます。

Artemisで疑似乱数を生成する場合はTRandomNumberEventStoreを使えば良いです。

まずは以下のようなyamlファイルを作成します。

ここではファイル名をNDimSphere.yamlとして、steeringの中に作成します。

Anchor:
  - &maxevt 1000000
  - &output output/sim/NDimSphere.root
Processor:
  - name: timer
    type: art::TTimerProcessor
##################################################################################################
  - name: MyTRandomNumberEventStore
    type: art::TRandomNumberEventStore
    parameter:
      Max: 1  # [Float_t] the maximum value
      MaxLoop: *maxevt  # [Int_t] the maximum number of loop
      Min: 0  # [Float_t] the minimum value
      OutputCollection: dim1  # [TString] output name of random values
      OutputTransparency: 0  # [Bool_t] Output is persistent if false (default)
      Verbose: 1  # [Int_t] verbose level (default 1 : non quiet)
##################################################################################################
  - name: MyTRandomNumberEventStore
    type: art::TRandomNumberEventStore
    parameter:
      Max: 1  # [Float_t] the maximum value
      MaxLoop: *maxevt  # [Int_t] the maximum number of loop
      Min: 0  # [Float_t] the minimum value
      OutputCollection: dim2  # [TString] output name of random values
      OutputTransparency: 0  # [Bool_t] Output is persistent if false (default)
      Verbose: 1  # [Int_t] verbose level (default 1 : non quiet)
##################################################################################################
  - name: MyTRandomNumberEventStore
    type: art::TRandomNumberEventStore
    parameter:
      Max: 1  # [Float_t] the maximum value
      MaxLoop: *maxevt  # [Int_t] the maximum number of loop
      Min: 0  # [Float_t] the minimum value
      OutputCollection: dim3  # [TString] output name of random values
      OutputTransparency: 0  # [Bool_t] Output is persistent if false (default)
      Verbose: 1  # [Int_t] verbose level (default 1 : non quiet)
##################################################################################################
  - name: MyTRandomNumberEventStore
    type: art::TRandomNumberEventStore
    parameter:
      Max: 1  # [Float_t] the maximum value
      MaxLoop: *maxevt  # [Int_t] the maximum number of loop
      Min: 0  # [Float_t] the minimum value
      OutputCollection: dim4  # [TString] output name of random values
      OutputTransparency: 0  # [Bool_t] Output is persistent if false (default)
      Verbose: 1  # [Int_t] verbose level (default 1 : non quiet)
##################################################################################################
  - name: MyTOutputTreeProcessor
    type: art::TOutputTreeProcessor
    parameter:
      FileName: *output  # [TString] The name of output file
      OutputTransparency: 0  # [Bool_t] Output is persistent if false (default)
      SplitLevel: 0  # [Int_t] Split level of tree defined in TTree (default is changed to be 0)
      TreeName: tree  # [TString] The name of output tree
      Verbose: 1  # [Int_t] verbose level (default 1 : non quiet)
##################################################################################################

次にartemisを起動します。

VNCを使う場合はDISPLAY=:`cat $vncdisplay` artemis -l $*を、使用しない場合はartemis -lを叩けば起動できます。

artemisが起動したら、add steering/NDimSphere.yaml、resを順番に実行し、解析を開始します。

TTimerProessorを入れているので解析が終わると以下のメッセージが出てきます。

Info in <art::TTimerProcessor::PostLoop>: real = 5.25, cpu = 5.21 sec, total 1000000 events, rate 191938.58 evts/sec

fcd 0でtreeのあるディレクトリに移動します。

正しく移動できていれば、brコマンドでyamlファイルで設定した名前でブランチができているはずです。

artemis [4] fcd 0
artemis [4] br
dim1                 art::TSimpleData    
dim2                 art::TSimpleData    
dim3                 art::TSimpleData    
dim4                 art::TSimpleData

モンテカルロ積分では領域内の被積分関数の期待を求める計算手法です。

N次元球の場合は、

$$ \mathrm{condition} = \left| \sum_{i=1}^{N} {r_i}^2 \right| < 1 $$

なのでtree->Draw("branch_name","condition","")で結果を描画することができます。

tree->Draw("dim1.fValue>>hypersphere1(100,0,1)","abs(pow(dim1.fValue,2))<1","colz")
tree->Draw("dim1.fValue:dim2.fValue>>hypersphere2(100,0,1,100,0,1)","abs(pow(dim1.fValue,2)+pow(dim2.fValue,2))<1","colz")
tree->Draw("dim1.fValue:dim2.fValue>>hypersphere3(100,0,1,100,0,1)","abs(pow(dim1.fValue,2)+pow(dim2.fValue,2)+pow(dim3.fValue,2))<1","colz")
tree->Draw("dim1.fValue:dim2.fValue>>hypersphere4(100,0,1,100,0,1)","abs(pow(dim1.fValue,2)+pow(dim2.fValue,2)+pow(dim3.fValue,2)+pow(dim4.fValue,2))<1","colz")
tree->Draw("dim1.fValue:dim2.fValue:dim3.fValue:dim4.fValue>>hypersphere4_4D","abs(pow(dim1.fValue,2)+pow(dim2.fValue,2)+pow(dim3.fValue,2)+pow(dim4.fValue,2))<1","colz")

$N=1$	$N=2$	$N=3$	$N=4$

本題の積分を計算していきます。

発生させた乱数の個数(loopを回した回数)と条件に入った個数の商が積分値です。

今回は区間が $[ 0,1 ]$なのでファクターとして $2^N$をかける必要があります。

analysisInfoにはloopを回した回数などの解析の情報が入っています。

今回使いたいのはloop数なのでanalysisInfo->GetAnalyzedEventNumber()を使います。

ヒント

analysisInfo->ls()でメソッドの一覧が見れます。

また条件に入った個数は、先ほど作成した図のエントリー数に対応していますので、histogram_name->GetEntries()でOKです。

よって以下のコマンドで計算できます。

hypersphere1->GetEntries()/analysisInfo->GetAnalyzedEventNumber()*2.
hypersphere2->GetEntries()/analysisInfo->GetAnalyzedEventNumber()*4.
hypersphere3->GetEntries()/analysisInfo->GetAnalyzedEventNumber()*8.
hypersphere4->GetEntries()/analysisInfo->GetAnalyzedEventNumber()*16.

これにて計算できました。

ちなみにちゃんと計算した時の値は以下のように計算できます。

Double_t r=1;
pow(r,1)*pow(TMath::Pi(),0)*2.
pow(r,2)*pow(TMath::Pi(),1)*1.
pow(r,3)*pow(TMath::Pi(),1)*(4./3.)
pow(r,4)*pow(TMath::Pi(),2)*(1./2.)

最後に比較をしてみます。

Double_t r=1;
abs(hypersphere1->GetEntries()/analysisInfo->GetAnalyzedEventNumber()*2.-pow(r,1)*pow(TMath::Pi(),0)*2.)
abs(hypersphere2->GetEntries()/analysisInfo->GetAnalyzedEventNumber()*4.-pow(r,2)*pow(TMath::Pi(),1)*1.)
abs(hypersphere3->GetEntries()/analysisInfo->GetAnalyzedEventNumber()*8.-pow(r,3)*pow(TMath::Pi(),1)*(4./3.))
abs(hypersphere4->GetEntries()/analysisInfo->GetAnalyzedEventNumber()*16.-pow(r,4)*pow(TMath::Pi(),2)*(1./2.))

返ってくる値を見ればわかりますが、3次元の時点で精度が悲しい感じになっています。

loop(yamlファイルの2行目maxevtの値)を変えて解析しなおしたり、一次元の図を眺めてみたりなどして原因を考えてみましょう。