入力
電磁場取り込み
User elemetns
結果処理
Python
GPT計算結果を処理する方法
GPTを使った計算結果の処理方法のメモです．
目次
基本
計算結果のデータ処理
GDF2A関係
特定の時刻/場所のデータの抽出
gdfファイルから，ある特定の時刻(time)や場所(position)の粒子のデータを取り出すことができず不便である．GDF2Aによりテキストに変換してから，GPTの特定の時刻や場所の計算結果を他のプログラムで処理するときに，そのことを思う．GDF2Aを使うと，全ての時刻と場所のデータがファイルに書かれ，とても巨大なファイルができあがるので，処理が大変である．
この問題を解決するために，必要な時刻あるいは場所の粒子のデータをファイルに書き出すためのPerlのプログラム(gdf2a_CutData.pl)を書いた．処理の流れは次のようにする．
- GDF2Aにより，gdfファイルをテキストに変換し，標準出力に出力する．
- GDF2Aの標準出力データをPerlのプログラム(gdf2a_CutData.pl)は読み込み，指定されたデータのみ標準出力に出力する．
- Perlのプログラムの標準出力は，リダイレクトによりファイルに保存する．
実際のコマンドは，パイプラインやリダイレクトで接続し，次のようにする．これは，Windows用のコマンド．
gdf2a -w 15 acc.gdf x y z G | perl gdf2a_CutData.pl position 0.8 0.9 > acc_part.dat
これで，positionが 0.8から 0.9 のデータのみが，テキストファイル(acc_part.dat)に書かれる．ここで使っているPerlのプログラム(gdf2a_CutData.pl)は，次のとおり．
計算結果の一部(時間，場所)を切り出すプログラム(gdf2a_CutData.pl)．
001 while(<STDIN>){ 002 @key_word=split(/\s+/); 003 chomp; 004 if($key_word[0]=~/$ARGV[0]/ && 005 $ARGV[1]<=$key_word[1] && $key_word[1]<=$ARGV[2]){ 006 007 print "\# ".$_."\n"; 008 009 while(<STDIN>){ 010 if($_=~/position/ || $_=~/time/){ 011 $read_data="TRUE"; 012 last; 013 }else{ 014 print $_; 015 } 016 } 017 018 } 019 020 if($read_data=~/TRUE/){ 021 $read_data=""; 022 redo; 023 } 024 }
あるエネルギー幅の粒子数
問題
加速器のビーム解析では，あるエネルギー幅に含まれる電流値を解析したい場合がある．ビームのエネルギー分布が下図のような場合，エネルギー幅ΔEのビーム電流は?—と言うような問題である．この問題では，図中のE0を変化させて，積分量の最大値を捜す必要がある．
|
| ビームのエネルギー分布．エネルギー幅ΔEの電流を計算したい． |
この問題は，特定のエネルギー範囲にある粒子数をカウントすることに置き換えることができる．図中のエネルギー範囲になる粒子数をカウントして，それぞれの電流値の対応する重み付けを乗算したものを積算する．すると，特定のエネルギー範囲にある電流値が分かる．後は，E0を変化させて，電流値の最大値を捜せばよい．
GPTには，このような問いに答えるコマンドは用意されていない．GPTのプログラムは非常に柔軟にできており，ユーザーがコマンドを追加し，この計算を行うことも可能である．しかし，GPTにこの機能を追加させることは思ったよりも面倒なので，ここでは専用の解析プログラムを示す．
プログラム
ソースプログラム
ソースプログラムはC++で書いている．私は，Linuxでコンパイル・実行を行った．普通のC++なので，Windowsでも同じように動作させることができるであろう．
メイン関数(main.cpp)．
001 #include "calEnergySpread.h" 002 #include <cstdlib> 003 004 int main(int argc, char *argv[]) 005 { 006 007 if(argc != 4){ 008 cerr << "Error! argument filename bandwidth precision" << endl; 009 exit(1); 010 } 011 012 cout << "filename:" << argv[1] << endl; 013 cout << "band width[KeV]:" << atof(argv[2])/1e3 014 << "\tprecision[KeV]:" << atof(argv[3])/1e3 << endl; 015 016 calEnergySpread *es = new calEnergySpread(argv[1]); 017 es->search_max_numpar(atof(argv[2]), atof(argv[3])); 018 019 cout << "E range[Mev]\t" 020 << es->peak_left/1e6 << "\t\t" 021 << es->peak_right/1e6 << "\t\t" 022 << "count:\t" << es->max_count << endl << endl << endl;; 023 024 return 0; 025 }
ヘッダーファイル(calEnergySpread.h)．
001 #ifndef CALENERGYSPREAD_H 002 #define CALENERGYSPREAD_H 003 004 #include <iostream> 005 #include <vector> 006 007 using namespace std; 008 009 class calEnergySpread{ 010 public: 011 calEnergySpread(char*); 012 ~calEnergySpread(); 013 int searchE_range(const double low, const double high); 014 void search_max_numpar(const double bw, const double step); 015 double *arrayE; 016 double min_E; 017 double max_E; 018 int max_count; 019 double peak_left; 020 double peak_right; 021 022 private: 023 vector<double> energy; 024 int numpar; 025 }; 026 027 #endif
ビーム電流値を計算するクラスの定義(calEnergySpread.cpp)．
001 #include "calEnergySpread.h" 002 #include <fstream> 003 #include <stdlib.h> 004 005 int compare(const void *a, const void *b); 006 007 //================================================================= 008 // member function calEnergySpread 009 //================================================================= 010 011 //================================================================= 012 // constructor 013 //================================================================= 014 calEnergySpread::calEnergySpread(char *filename) 015 { 016 ifstream datain(filename); 017 018 if(!datain){ 019 cerr << "can not open a file. file-name :" << filename << endl; 020 exit(1); 021 } 022 023 //--- 二行のメッセージを読み飛ばす --- 024 string temp; 025 getline(datain, temp); 026 getline(datain, temp); 027 028 unsigned int in=0; 029 double x, y, z, G; 030 031 while(!((datain >> x >> y >> z >> G).eof())){ 032 if(in>=energy.size()) energy.resize(energy.size()+100); // サイズが不足すると拡張 033 energy[in] = 511e3*(G-1); 034 in++; 035 } 036 037 numpar = in; 038 datain.close(); 039 energy.resize(in); 040 041 cout << "numper = " << numpar << "\t size="<<energy.size() << endl; 042 043 arrayE = new double [energy.size()]; 044 045 for(unsigned int i=0; i<energy.size(); i++){ 046 arrayE[i]=energy[i]; 047 } 048 049 qsort(arrayE, energy.size(), sizeof(double), compare); 050 051 cout << "lowest = " << arrayE[0] << endl; 052 cout << "highest = " << arrayE[energy.size()-1] << endl; 053 054 min_E = arrayE[0]; 055 max_E = arrayE[energy.size()-1]; 056 057 } 058 059 060 //================================================================= 061 // constructor 062 //================================================================= 063 calEnergySpread::~calEnergySpread() 064 { 065 delete arrayE; 066 } 067 068 069 //================================================================= 070 // search energy range 071 //================================================================= 072 int calEnergySpread::searchE_range(const double low, const double high) 073 { 074 static int left, right, middle; 075 static int lower, higher; 076 077 //------- lower limitを探す -- 078 left=0; 079 right=energy.size()-1; 080 081 middle=(left+right)/2; 082 083 while(right-left>1){ 084 if(arrayE[middle] <= low){ 085 left=middle; 086 }else{ 087 right=middle; 088 } 089 middle=(left+right)/2; 090 } 091 092 lower=right; 093 094 //------- higher limitを探す -- 095 left=0; 096 right=energy.size()-1; 097 098 middle=(left+right)/2; 099 100 while(right-left>1){ 101 if(arrayE[middle] < high){ 102 left=middle; 103 }else{ 104 right=middle; 105 } 106 middle=(left+right)/2; 107 } 108 109 higher=left; 110 111 return higher-lower+1; 112 113 114 } 115 116 //================================================================= 117 // search energy range 118 //================================================================= 119 void calEnergySpread::search_max_numpar(const double bw, const double step) 120 { 121 static unsigned int numarray; 122 static unsigned int max_i; 123 124 numarray=(int)((max_E-min_E)/step-bw/step)+1; 125 126 int *count = new int [numarray]; 127 128 for(unsigned int i=0; i<numarray; i++){ 129 count[i]=searchE_range(min_E+i*step, min_E+i*step+bw); 130 } 131 132 max_count = 0; 133 for(unsigned int i=0; i<numarray; i++){ 134 if(max_count < count[i]){ 135 max_count = count[i]; 136 max_i=i; 137 } 138 } 139 140 peak_left = min_E + max_i*step; 141 peak_right = min_E + max_i*step+bw; 142 143 } 144 145 //================================================================= 146 // compare function for qsort 147 //================================================================= 148 int compare(const void *a, const void *b) 149 { 150 151 if(*(double *)a < *(double *)b){ 152 return -1; 153 }else{ 154 return 1; 155 } 156 157 }
コンパイル方法
このMakefileを使って，makeで実行ファイルができる．
使い方
このプログラムは，GPTのある特定の時刻，あるいは場所でのあるエネルギー範囲に含まれている粒子数を数える．特定の時刻/場所の粒子のデータは，先に述べた特定の時刻/場所のデータの抽出を使う．具体的な計算手順は，以下の通り．
- GPTによりビームの軌道を計算する．このとき，時刻，あるいは場所でビーム情報を出力(gdf)する必要がある．
- 粒子数をカウントしたい特定の時刻/場所のビーム情報のファイルを作る(特定の時刻/場所のデータの抽出を参照)．このときのファイルは，ひとつの時刻，あるいは場所でなくてはならない．一つのファイルに，複数の時刻/場所の情報が書き込まれると，エラーあるいは誤った計算結果が出力される．
- 作成したプログラムに，必要なコマンドライン引数を指定して実行させさせる．
calES ファイル名 エネルギー幅[eV] エネルギー幅の精度[eV]
calESが実行ファイル名で，引き続き解析するファイル名，粒子数をカウントするエネルギー幅ΔE，E0を変化させるときの刻み幅を指定する．解析するデータが多くなると，シェルスクリプトを書いた方が良い(シェルスクリプトの例)．そして，計算結果をリダイレクトでファイルに保存する．シェルスクリプト > result.dat