2014-02-18

ROOT Treeからヒストグラム

研究 ROOT

t1->Draw("a[0]:a[1]");

から、

t1->Draw("a[0]:a[1]","","COLZ");

ってやると２次元ヒストがヒートマップみたいになる

2014-02-06

ROOTでコンパイル

研究 ROOT

ROOTでコンパイルするためのMakefile

TARGETS=Hello

ROOTFLAGS = $(shell root-config --cflags)
ROOTLIBS  = $(shell root-config --libs)

CXXFLAGS  = -Wall -O2 $(ROOTFLAGS)
CXXLIBS   = $(ROOTLIBS)

all: $(TARGETS)

Hello: Hello.o
        g++ -o $@ Hello.o $(CXXLIBS)

.cc.o:
        g++ -c $(CXXFLAGS) $<

2014-01-23

vector便利

プログラミング C++

配列って最初定義するときに要素数が必要だけど、vectorならいらない。

基本的な使い方は、イテレータと合わせて

#include <vector>
#include <iostream>

using namespace std;

int main(){
    vector<int> array;
    int i;

    for(i=0;i<10;i++){
        array.push_back(i);
    }

    for(i=0;i<10;i++){
        cout << array[i] << endl;
    }

    vector<int>::iterator it = array.begin();

    while(it != array.end()){
        cout << *it << endl;
        ++it;
    }
}

2014-01-23

gslでできること一覧

gsl 物理学研究

gslでできること一覧(マニュアルより抜粋)

複素数多項式の求根法特殊関数置換組み合わせ整列
線形代数 CBLAS ライブラリ高速フーリエ変換乱数数値積分
乱数分布ヒストグラム統計モンテカルロ積分微分方程式
シミュレーティド・アニーリングベクトルと行列 BLAS の利用
固有値問題疑似乱数列 N項組数値微分補間級数展開
チェビシェフ近似求根法離散ハンケル変換最小二乗法最適化
ウェーブレット IEEE 浮動小数点表現物理定数

やりたいことがあれば、マニュアルで関数を見てやるかんじ

2014-01-23

gslの使い方

物理学研究 gsl

gslとは、GNU Scientific Library。いろいろめんどくさい計算すぐできる。

インストールしただけではただのライブラリ。

例：ベッセル関数 \(J_0(x)\) の、\(x=5\) での値を計算して出力する

#include <stdio.h>
#include <gsl/gsl_sf_bessel.h>

int main (void){
          double x = 5.0;
          double y = gsl_sf_bessel_J0(x);
          printf("J0(%g) = %.18e\n", x, y);
          return 0;
}

として、コンパイルは

$ gcc $(gsl-config --cflags) example.c $(gsl-config --libs)

とする。

ちなみに、

$ echo $(gsl-config --cflags)
-I/usr/local/include
$ echo $(gsl-config --libs)
-L/usr/local/lib -lgsl -lgslcblas -lm

で、できたa.outを実行すると、

$ ./a.out
J0(5) = -1.775967713143382642e-01

2014-01-22

Geant4でのデータ収集の仕方(2)

研究 Geant4 物理学

一番簡単なやり方、

ユーザフック

exampleB1などで使われる。B1では線量計算をStep、Event、Run actionを自作して求める。

1. 構造物を作り、stepping actionをするvolumeを指定
2. Stepping actionクラス内で、今のStepがスコアすべき物体中にあるかどうかを調べて、その場合そのStepでのエネルギー付与を求めて積算
3. Event actionクラスで各Eventの後でエネルギー付与の積算
4. Run actionクラスで全Eventの積算を求めて、スコアした物体の質量と合わせて線量[J/kg]を計算し、出力　

具体的に、
1. 構造物を作り、stepping actionをするvolumeを指定

new G4PVPlacement(0,                       //no rotation
                  pos2,                    //at position
                  logicShape2,             //its logical volume
                  "Shape2",                //its name
                  logicEnv,                //its mother  volume
                  false,                   //no boolean operation
                  0,                       //copy number
                  checkOverlaps);          //overlaps checking

// Set scoring volume to stepping action
// (where we will account energy deposit)
//
B1SteppingAction* steppingAction = B1SteppingAction::Instance();
////steppingAction->SetVolume(logicShape1);
steppingAction->SetVolume(logicShape2);

2. Stepping actionクラス内で、今のStepがスコアすべき物体中にあるかどうかを調べて、その場合そのStepでのエネルギー付与を求めて積算

void B1SteppingAction::UserSteppingAction(const G4Step* step)
{
  // get volume of the current step
  G4LogicalVolume* volume
    = step->GetPreStepPoint()->GetTouchableHandle()
      ->GetVolume()->GetLogicalVolume();

  // check if we are in scoring volume
  if (volume != fVolume ) return;

  // collect energy and track length step by step
  G4double edep = step->GetTotalEnergyDeposit();
  fEnergy += edep;
}

3. Event actionクラスで各Eventの後でエネルギー付与の積算

void B1EventAction::BeginOfEventAction(const G4Event* event)
{
  G4int eventNb = event->GetEventID();
  if (eventNb%fPrintModulo == 0) {
    G4cout << "\n---> Begin of event: " << eventNb << G4endl;
  }

  // Reset accounted energy in stepping action
  B1SteppingAction::Instance()->Reset();
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void B1EventAction::EndOfEventAction(const G4Event* /*event*/)
{
  // accumulate statistics
  G4double energy = B1SteppingAction::Instance()->GetEnergy();
  fEnergySum  += energy;
  fEnergy2Sum += energy*energy;
}

4. Run actionクラスで全Eventの積算を求めて、スコアした物体の質量と合わせて線量[J/kg]を計算し、出力

// Compute dose
//
G4double energySum  = B1EventAction::Instance()->GetEnergySum();
G4double energy2Sum = B1EventAction::Instance()->GetEnergy2Sum();
G4double rms = energy2Sum - energySum*energySum/nofEvents;
if (rms > 0.) rms = std::sqrt(rms); else rms = 0.;

G4double mass = B1SteppingAction::Instance()->GetVolume()->GetMass();
G4double dose = energySum/mass;
G4double rmsDose = rms/mass;

…

// Print
//
G4cout
   << "\n--------------------End of Run------------------------------\n"
   << " The run consists of " << nofEvents << " "<< particleName << " of "
   <<   G4BestUnit(particleEnergy,"Energy")
   << "\n Dose in scoring volume "
   << B1SteppingAction::Instance()->GetVolume()->GetName() << " : "
   << G4BestUnit(dose,"Dose")
   << " +- "                   << G4BestUnit(rmsDose,"Dose")
   << "\n------------------------------------------------------------\n"
   << G4endl;

2014-01-22

Geant4でのデータ収集の仕方(1)

Geant4 C++ 物理学研究

Geant4でデータを取る(スコアリング)ために用意されているもの

Hit

検出器のsensitiveな領域内での、Trackの物理相互作用のスナップショット。G4Stepオブジェクトに属する以下のような物理量を集める

そのStepの位置と時刻
Trackの運動量、エネルギー
そのStepでのエネルギー付与
ジオメトリ情報

Hits Collection

G4Stepに属するもの

Sensitive Detector

G4Stepオブジェクトを入力として、StepのHit情報の集合を作る。
検出器の中で指定するには、VolumeにSensitive Detectorへのポインタを持たせる

例(expample B4c)

B4cCalorimeterSD* absoSD
    = new B4cCalorimeterSD("AbsorberSD", "AbsorberHitsCollection", nofLayers); 

G4SDManager::GetSDMpointer()->AddNewDetector(absoSD );
absorberLV->SetSensitiveDetector(absoSD);

（旧）研究メモ

kennkyuumemo

ROOT Treeからヒストグラム

ROOTでコンパイル

vector便利

gslでできること一覧

gslの使い方

Geant4でのデータ収集の仕方(2)

ユーザフック

Geant4でのデータ収集の仕方(1)

Hit

Hits Collection

Sensitive Detector