物理の駅 by onsanai

Physics-station 研究で日々感じたことを忘れないための備忘録

Anaconda updateでのエラー

conda update -all

で最後に

Preparing transaction: done
Verifying transaction: done
Executing transaction: | DEBUG menuinst_win32:__init__(196): Menu: name: 'Anaconda${PY_VER} ${PLATFORM}', prefix: 'C:\Users\%username%\Anaconda3', env_name: 'None', mode: 'user', used_mode: 'user'
DEBUG menuinst_win32:create(320): Shortcut cmd is C:\Users\%username%\Anaconda3\python.exe, args are ['C:\\Users\\%username%\\Anaconda3\\cwp.py', 'C:\\Users\\%username%\\Anaconda3', 'C:\\Users\\%username%\\Anaconda3\\python.exe', 'C:\\Users\\%username%\\Anaconda3\\Scripts\\jupyter-notebook-script.py', '%USERPROFILE%']
/ DEBUG menuinst_win32:__init__(196): Menu: name: 'Anaconda${PY_VER} ${PLATFORM}', prefix: 'C:\Users\%username%\Anaconda3', env_name: 'None', mode: 'user', used_mode: 'user'
DEBUG menuinst_win32:create(320): Shortcut cmd is C:\Users\%username%\Anaconda3\python.exe, args are ['C:\\Users\\%username%\\Anaconda3\\cwp.py', 'C:\\Users\\%username%\\Anaconda3', 'C:\\Users\\%username%\\Anaconda3\\python.exe', 'C:\\Users\\%username%\\Anaconda3\\Scripts\\jupyter-notebook-script.py', '%USERPROFILE%']
done

というエラーが出るのはなぜだろう

Slackのあるチャンネルの投稿を全て消すpythonコード

import time
from slackclient import SlackClient
slack_client = SlackClient('****-************-************-************-********************************')
target_channel ="channel"

def list_channels():
    channels_call = slack_client.api_call("channels.list")
    if channels_call.get('ok'):
        return channels_call['channels']
    return None

def channel_info(channel_id):
    channel_info = slack_client.api_call("channels.info", channel=channel_id)
    if channel_info:
        return channel_info['channel']
    return None

if __name__ == '__main__':
    channels = list_channels()
    if channels:
        for channel in channels:
            print(channel['name'])
            if channel['name'] != target_channel:
                continue            
            while True:
                detailed_info = channel_info(channel['id'])
                if not detailed_info['latest']:
                    break                
                message_ts = detailed_info['latest']['ts']
                slack_client.api_call("chat.delete", channel=channel['id'], ts=message_ts)
                time.sleep(1)
    else:
        print("Unable to authenticate.")

削除は出来るが、Total messagesの数は減らない。

OKIのプリンターでトナーの使う量を減らす(節約する)

OKI C811の場合

f:id:onsanai:20180613144845p:plain:w400

ログイン -> 管理者設定 -> 印刷設定 -> トナーセーブ -> トナーセーブ量 -> やや多い

選択肢は オフ・少ない・やや多い・多いの4つある

OKI B841の場合

f:id:onsanai:20180613145334p:plain:w400

ログイン -> プリンタ -> 印刷メニュー -> 印刷品質 -> トナーセーブモード -> オン

選択肢はユウコウ・ムコウの2通り

ちなみに初期パスワードは「aaaaaa」

WindowsでVPNを使うために必須の機能 RAS 接続マネージャー管理キット (CMAK)

今まで使えていたVPNが突然使えなくなった。

ユーザー は終了した  という接続をダイヤルしました。終了時に戻された理由コードは 829 です。

が出てVPNが接続できなくなったが、原因は Windowsの機能の有効化または無効化で

RAS 接続マネージャー管理キット (CMAK)

をオフにしたのが原因だった。オンにして再起動でVPNが使えるようになった。

ツイッターで検索したときにキーボードショートカットで「j、k」でツイートを選択したときにツイート内容が隠れるバグ

ツイッター公式には、ホーム画面や検索画面でキーボードショートカットを使うことが出来る。どういうショートカットか有効かは「?」を画面で押してみると良い。

f:id:onsanai:20180610014929p:plain:w500

ホーム画面では、このショートカットで次のツイートを表示する「j」前のツイートを表示する「k」はよく働いている。しかし、検索画面でこの機能を使うと、どういう内容で検索したかというバーが邪魔をし、選択したツイートが隠れてしまうと言うバグが存在する。この1年ほどこのバグが修正されていないのだけれど、だれか解決方法知っている人はいませんか?

Windows10でシステムの冷却ポリシーや最大のプロセッサの状態を表示させる

ascii.jp

この内容はWindows10でも有効である。

Surface Pro 3 のOSのアップデートをすると、以前有効にした電源の詳細設定の項目がなくなってしまう。上記の記事を参考に、以下の項目の Attributes キーを作成し、 DWORD2 を追加した。

  • 最大のプロセッサの状態 bc5038f7-23e0-4960-96da-33abaf5935ec
  • システムの冷却ポリシー 94D3A615-A899-4AC5-AE2B-E4D8F634367F

冷却ポリシーは、ざっくり言えばファンの挙動を変えてくれる。

f:id:onsanai:20180528100303p:plain

表示されましたね。

Windows 10 で 3Dオブジェクトフォルダを削除する

アップグレードのたびにフォルダが作成されるのがウザいので、レジストリから削除するスクリプトを作った。 delete_3Dobj.reg とでもして管理者権限で実行されたし。

Windows Registry Editor Version 5.00

[-HKEY_LOCAL_MACHINE\SOFTWARE\Microsoft\Windows\CurrentVersion\Explorer\MyComputer\NameSpace\{0DB7E03F-FC29-4DC6-9020-FF41B59E513A}]
[-HKEY_LOCAL_MACHINE\SOFTWARE\Microsoft\Windows\CurrentVersion\Explorer\MyComputer\NameSpace\{3dfdf296-dbec-4fb4-81d1-6a3438bcf4de}]
[-HKEY_LOCAL_MACHINE\SOFTWARE\Microsoft\Windows\CurrentVersion\Explorer\MyComputer\NameSpace\{24ad3ad4-a569-4530-98e1-ab02f9417aa8}]
[-HKEY_LOCAL_MACHINE\SOFTWARE\Microsoft\Windows\CurrentVersion\Explorer\MyComputer\NameSpace\{f86fa3ab-70d2-4fc7-9c99-fcbf05467f3a}]
[-HKEY_LOCAL_MACHINE\SOFTWARE\Microsoft\Windows\CurrentVersion\Explorer\MyComputer\NameSpace\{d3162b92-9365-467a-956b-92703aca08af}]

デスクトップ {B4BFCC3A-DB2C-424C-B029-7FE99A87C641} ダウンロード {088e3905-0323-4b02-9826-5d99428e115f}

は残したいので入れず。

全部消す場合は以下のように。上から順に3Dオブジェクト ドキュメント 音楽 ピクチャー ビデオ デスクトップ ダウンロードです。

Windows Registry Editor Version 5.00

[-HKEY_LOCAL_MACHINE\SOFTWARE\Microsoft\Windows\CurrentVersion\Explorer\MyComputer\NameSpace\{0DB7E03F-FC29-4DC6-9020-FF41B59E513A}]
[-HKEY_LOCAL_MACHINE\SOFTWARE\Microsoft\Windows\CurrentVersion\Explorer\MyComputer\NameSpace\{3dfdf296-dbec-4fb4-81d1-6a3438bcf4de}]
[-HKEY_LOCAL_MACHINE\SOFTWARE\Microsoft\Windows\CurrentVersion\Explorer\MyComputer\NameSpace\{24ad3ad4-a569-4530-98e1-ab02f9417aa8}]
[-HKEY_LOCAL_MACHINE\SOFTWARE\Microsoft\Windows\CurrentVersion\Explorer\MyComputer\NameSpace\{f86fa3ab-70d2-4fc7-9c99-fcbf05467f3a}]
[-HKEY_LOCAL_MACHINE\SOFTWARE\Microsoft\Windows\CurrentVersion\Explorer\MyComputer\NameSpace\{d3162b92-9365-467a-956b-92703aca08af}]
[-HKEY_LOCAL_MACHINE\SOFTWARE\Microsoft\Windows\CurrentVersion\Explorer\MyComputer\NameSpace\{B4BFCC3A-DB2C-424C-B029-7FE99A87C641}]
[-HKEY_LOCAL_MACHINE\SOFTWARE\Microsoft\Windows\CurrentVersion\Explorer\MyComputer\NameSpace\{088e3905-0323-4b02-9826-5d99428e115f}]

参照

blog.yuizi.com

Windows 10 Fall Creators Updateにて現れた 3Dオブジェクトフォルダを削除する – 友蔵保管庫

Geant4で原子核乾板 (10の12乗チャンネル) を実装した

ミクロンの空間分解能を持ちながら、ミリメートルからメートルのサイズを持つ放射線検出器は原子核乾板以外に存在しない。通常、チャンネル数はある特定の数になるが、原子核乾板の場合はそういう概念はないに近い。1平方センチあたり10の14乗を程度チャンネル数が存在するとみなしてもよい。

高エネルギー粒子をシミュレーションするときは、乾板の表面に入ってから出ていくまでを直線と仮定して差し支えないが、低エネルギー粒子の飛跡をシミュレーションで再現しようとするときは、PhysicsListが低エネルギー側でどの程度現実の物理を反映しているかどうかは別にして、それなりに細かいStepで飛跡の位置情報を出力したいと誰しも思うだろう。

Geant4でこのような検出器を実装するには、 G4PVReplica という機能を使う。

G4PVReplicaG4VPhysicalVolume から継承されているクラスで、先に作った箱(例えば G4Box)をある方向に分割する機能である。X軸、Y軸、Z軸の3回分割を繰り返せば、細かい要素に箱を区切ることができる。シミュレーションしたい内容に応じて、細かさは調整すればよい。 コメントで指摘があったとおり、最小要素(Elmt)を四角柱(Stick)に詰め、それを平面に伸びる四角柱(Srfc)に詰め、さらに箱(Box)に詰めている、という考え方の方が分かりやすい。

次の例では、簡単のために一辺100ミリメートルの立方体を、一辺10ミクロンの立方体で区切る。計10の12乗のチャンネル数がGeant4にできていることになる。

G4double leng_EmulsionBox = 100.0*mm;     // X-full-length of emulsion box 
G4double leng_EmulsionElmt = 0.01*mm;     // X-full-length of emulsion element
G4int nDiv = std::round(leng_EmulsionBox / leng_EmulsionElmt);

//Emulsionを格納するBOX
auto* solid_EmulsionBox =
    new G4Box("Solid_EmulsionBox", leng_EmulsionBox / 2.0, leng_EmulsionBox / 2.0, leng_EmulsionBox / 2.0);
auto* materi_EmulsionBox = materi_Man->FindOrBuildMaterial("G4_AIR");
auto* logVol_EmulsionBox =
    new G4LogicalVolume(solid_EmulsionBox, materi_EmulsionBox, "LogVol_EmulsionBox", 0, 0, 0);

//XYに伸びる四角柱(平面)
auto* solid_EmulsionSrfc =
    new G4Box("Solid_EmulsionSrfc", leng_EmulsionBox / 2.0, leng_EmulsionBox / 2.0, leng_EmulsionElmt / 2.0);
auto* materi_EmulsionSrfc = materi_Man->FindOrBuildMaterial("G4_AIR");
auto* logVol_EmulsionSrfc =
    new G4LogicalVolume(solid_EmulsionSrfc, materi_EmulsionSrfc, "LogVol_EmulsionSrfc", 0, 0, 0);
new G4PVReplica("PhysVol_EmulsionSrfc", logVol_EmulsionSrfc, logVol_EmulsionBox, kZAxis,
    nDiv, leng_EmulsionElmt);

//X方向に伸びる四角柱(棒)
auto* solid_EmulsionStick =
    new G4Box("Solid_EmulsionStick", leng_EmulsionBox / 2.0, leng_EmulsionElmt / 2.0, leng_EmulsionElmt / 2.0);
auto* materi_EmulsionStick = materi_Man->FindOrBuildMaterial("G4_AIR");
auto* logVol_EmulsionStick =
    new G4LogicalVolume(solid_EmulsionStick, materi_EmulsionStick, "LogVol_EmulsionStick", 0, 0, 0);
new G4PVReplica("PhysVol_EmulsionStick", logVol_EmulsionStick, logVol_EmulsionSrfc, kYAxis,
    nDiv, leng_EmulsionElmt);

//最小単位 (Element)
auto* solid_EmulsionElmt =
    new G4Box("Solid_EmulsionElmt", leng_EmulsionElmt / 2.0, leng_EmulsionElmt / 2.0, leng_EmulsionElmt / 2.0);
auto* materi_EmulsionElmt = materi_Man->FindOrBuildMaterial("EmulsionE373");
auto* logVol_EmulsionElmt =
    new G4LogicalVolume(solid_EmulsionElmt, materi_EmulsionElmt, "LogVol_EmulsionElmt", 0, 0, 0);
new G4PVReplica("PhysVol_EmulsionElmt", logVol_EmulsionElmt, logVol_EmulsionStick, kXAxis,
    nDiv, leng_EmulsionElmt);

auto threeVect_LogV = G4ThreeVector(0, 0, 0);
auto trans3D_LogV = G4Transform3D(G4RotationMatrix(), threeVect_LogV);
new G4PVPlacement(trans3D_LogV, "Emulsion", logVol_EmulsionBox, physVol_World, false, 0);

こういう使い方は、Geant4の開発者が想定する使用方法ではないと思うが、大丈夫ですか? (エラい人へ)

参考文献

https://wiki.kek.jp/display/geant4/Geant4+Japanese+Tutorial+for+Detector+Simulation+2017

原子核乾板の元素組成と密度

更新履歴を含めた情報は、Gitlabで管理しています

https://gitlab.com/yoshimoto/emulsiondensity/blob/master/ja.md

原子核乾板の密度と組成について

原子核乾板 (Nuclear Emulsion) の組成や密度を実測することは難しい。 軽元素から重元素までの測定レンジを持つ組成の測定方法は存在しないため、同じ条件で全ての元素組成を測定することはできない。 また、環境の湿度が変化すると乾板に含まれる水分量も変化し、酸素や水素の組成が変化する。 製造工程から推定する方法については、原料の組成は測定しやすいものの、原子核乳剤を製造する工程に水洗工程があり、この工程で何がどれだけ流出したかの推定は難しい。

よくある測定方法

軽元素である水素、炭素、窒素、硫黄、酸素の質量比(通常の装置では水素、炭素、窒素)は、元素分析装置による測定が可能である1。硫黄についてはその存在比が少ないため、元素分析装置での測定は一般には不可能である。 重元素については、銀、臭素ヨウ素の比はSEM-EDXによる測定が可能である[^8]。 定着(FIX)前後の厚みと質量測定により、AgBrIとバインダーの比及び密度を測定することが可能である[^8]。

アルファ線を用いた密度測定方法

アルファ線の飛程からエネルギーとレンジ(飛程)の関係式を導き出す方法が最も精度の良い密度測定方法である2

アルファ線の飛程は約1μmの標準偏差を持つ分布となり、100イベント使えば飛程につく標準誤差は0.1μmとなる。 密度3.5g/ccのPo-212からのアルファ線の飛程は約50μmなので、飛程につく誤差は0.2%となる。この誤差は、密度の誤差に換算すると0.3%相当になる。

この手法で求められる密度は、アルファ線のエネルギーとレンジ(飛程)の関係式を媒介する、換算密度とも呼ばれるものであり、絶対値としての信頼性ない。

原子核乾板の組成および密度

以下、論文等で紹介されている各原子核乾板の組成を紹介する。

Emulsion Ilford G-5

  • Density 3.907g/cc 3
H C N O S Ag Br I
Mass % ^21 1.4 6.9 1.7 6.9 0.2 47.4 34.8 0.6

50% R.H. 。

[^13]の第3表のIlford G-5 の Nが0.67g/ccとなっているのは0.067g/ccの誤記と思われる。4

Photographic Emulsion (ICRU-215) = Standard emulsion

  • Density 3.815g/cc 5
H C N O S Ag Br I
Mass % [^3] 1.4 7.2 1.9 6.6 0.2 47.4 34.9 0.3

Fuji ET-7B ET-7C/D

  • 密度: 3.73 g/cc (Fuji ET-7B)
  • 密度: 3.60 g/cc (Fuji ET-7C/D)
  • 密度: 3.63g/cc (ET-7D) [^11]
  • 相対湿度: 68%

結晶サイズ

  • ET-7B 240±78 nm 6
  • ET-7C 260 nm ^9
  • ET-7D 176 nm ^9
H C N O S Ag Br I
ET-7D Mass g/cc^9 0.05 0.31 0.1 0.226 0.0067 1.675 1.22 0.034
ET-7D Mass %^9 1.38 8.56 2.76 6.24 0.18 46.25 33.69 0.94
ET-7C/7D Mass %1 1.5 9.3 3.1 6.8 0.2 45.4 33.4 0.3
Binder AgBrI
ET-7D Volume %[^11] 55.0 45.0

Fuji OPERA film

  • 密度: 2.71g/cc 7
  • 密度: 2.84g/cc 乳剤層のみ [^11]
  • 密度: 2.77g/cc 乳剤層と保護層 [^11]

OPERAの論文(シミュレーションも多分)は全ては密度2.71g/cc と表記されている。桑原氏の値となぜズレたのかは不明。

結晶サイズ

  • 200±16 nm [^7]

組成

以下の組成は乳剤層と保護層である。

H C N O S Si Na Sr Ba Ag Br I
Mass % ^6 2.4 13.0 4.81 12.43 0.09 0.08 0.08 0.02 0.01 38.34 27.86 0.81
Binder AgBrI
Mass % 33.0 67.0
Volume %^7 69.0 31.0

組成は原料からの計算値である。水洗工程でバインダーが水溶液側に抜ける割合に大胆な仮定が入っている[^11]。

OPERAタイプ GIF処方[^11]

  • 結晶サイズ 200 nm
  • 密度 3.53g/cc
H C N O S Ag Br I
GIF g/cc[^11] 0.05 0.326 0.11 0.23 - 1.6 1.166 0.033
GIF Mass %[^11] 1.42 9.27 3.13 6.54 - 45.52 33.17 0.94

Nagoya NIT [^8]

  • 乳剤層の密度: 3.44 g/cc
  • AgBrの密度: 6.473g/cc
  • ゼラチンの密度: 1.32g/cc
  • PVAの密度: 1.19g/cc

結晶サイズ

  • NIT 44.2±6.8 nm
  • UNIT 24.8±4.3 nm

組成

H C N O Ag Br I
Mass % 1.6 10.1 2.7 7.4 44.5 31.8 1.9
Atomic % 41.1 21.4 4.9 11.7 10.5 10.1 0.4
Binder AgBrI
Mass % 21.9 78.2
Density g/cc 1.29 6.473

AgBrIの分量は原料からの計算値、HとCとNは実測値、Oはそこからの推定値である。

文責: 吉本

間違いがあればメール or コメントをください

参考文献


  1. 浅田 吉本ら PTEP論文 https://doi.org/10.1093/ptep/ptx076 (2017)

  2. Geant4で伝統的に使われてきたFujiFilm ET-7C/7D emulsion for KEK-PS E373 の組成 元ソースは歳藤氏による

  3. 古関 靖夫 原子核乾板 http://hdl.handle.net/2261/30908 (1959)

  4. 吉田ら NIM論文 https://doi.org/10.1016/j.nima.2016.11.044 (2016)

  5. SRIM ソフトウェア http://www.srim.org/

  6. 桑原謙一氏・吉本 private communication

  7. Photographic Emulsion (ICRU-215) https://physics.nist.gov/cgi-bin/Star/compos.pl?refer=ap&matno=215

  8. 中村琢ら NIM論文 https://dx.doi.org/10.1016/j.nima.2005.08.109 (2006)

  9. L. Patrizii ら JINST論文 https://dx.doi.org/10.1088/1748-0221/3/07/P07002 (2008)

直交座標から球面座標系(極座標)に変換するときの誤差伝搬

Error propagation when converting from rectangular coordinate system to spherical coordinate system

飛跡のベクトルが次のように得られたとする

Assume that a vector of trajectory is obtained as follows

 \displaystyle

V_x = x \pm \delta x \\
V_y = y \pm \delta y \\
V_z = z \pm \delta z \\

Range、Theta、Phiは次のようになる

Range, Theta, Phi are as follows

 \displaystyle
Range = r = \sqrt{x^2+y^2+z^2}\\
Theta = \theta =  acos(z/r)\\
Phi = \phi = atan(y/x)\\

誤差伝搬を計算すると誤差は次のようになる

By calculating the error propagation, the error becomes as follows

 \displaystyle
\delta r = \frac{ \sqrt{  (x\,\delta x)^2+ (y\,\delta y)^2+ (z\,\delta z)^2 }}{r} \\

\delta \theta = \frac { \sqrt{ (x^2 \, \delta x^2+y^2 \, \delta y^2)z^2      + (x^2+y^2)^2\delta z^2    }    }                {r^2\sqrt{x^2+y^2}}\\

\delta \phi =        \frac{\sqrt{(y \, \delta x)^2+(x \, \delta y)^2}}       {x^2+y^2}

検索用

参照

Derivative Calculator with Steps

 \displaystyle
x = r \sin\theta \cos\phi \\
y = r \sin\theta \sin\phi\\
z = r \cos\theta\\

Rangeを1に規格化するとき、

 \displaystyle
\delta x_{norm} = \sqrt{ (\cos\theta \cos\phi \cdot \delta \theta)^2+(\sin\theta\sin\phi\cdot\delta \phi)^2 }\\
\delta y_{norm} = \sqrt{ (\cos\theta \sin\phi \cdot\delta \theta)^2+(\sin \theta\cos \phi \cdot\delta \phi)^2}\\
\delta z_{norm} = \sqrt{(\sin\theta \cdot\delta \phi)^2}\\

Gitでサブモジュールの追加

新しいサブモジュールを追加。 root_macros フォルダに リポジトリ root_macros.gitを追加する。

git submodule add https://gitlab.com/yoshimoto/root_macros.git root_macros

リモートが更新されたら取り込む

git submodule foreach git pull origin master

[Git] git submodule は癖がすごいとの噂だったが素直につきあっていけそうという話 | deadwood git submoduleしてるリポジトリをリモートの最新に更新する - Qiita

Pattern matching sample パターンマッチングサンプル

(ja) パターンマッチング可能な飛跡ファイル

https://1drv.ms/f/s!Ap9xAxIuzM0xlLxu6slGIjkY3gSjNg

Q1. どういう手段でもいいので、 beam_4372-2_u.txt の飛跡と beam_4372-2_d.txt の飛跡の位置ずれを計算しよう。

Q2. 計算方法について、互いに紹介しよう。

Q3. これは、位置ずれのピークを求めるソースコードである。 std::map はどういう働きをするのか調べよう。

Q4. 前回配布したファイルはパターンマッチングが成功しない。成功しないと言うためにはどうすればよいか考えよう。

https://1drv.ms/t/s!Ap9xAxIuzM0xlLlaOr2g5qaKpYqF2w

https://1drv.ms/t/s!Ap9xAxIuzM0xlLlbUgqNsGvbe4Eq6w

(en)

Trajectory file for pattern-matching

https://1drv.ms/f/s!Ap9xAxIuzM0xlLxu6slGIjkY3gSjNg

Q1. Calculate the position shift between the trajectory of beam_4372-2_u.txt and the trajectory ofbeam_4372-2_d.txt with any way.

Q2. Introduce each other about your calculation methods.

Q3. This is a source code to search peak position shift . Examine what std::map does.

  • Position shift of all trajectories are packed in a histogram for each bin_width .
  • Create a blurred histogram with surrounding n_bin_bg x n_bin_bg .
  • Subtract the histogram from the blurred histogram
  • Find the peak with the the maximum content

Q4. Pattern matching does not succeed in files distributed last time. Let's see what we should say to show that the pattern matching does not succeed.

https://1drv.ms/t/s!Ap9xAxIuzM0xlLlaOr2g5qaKpYqF2w

https://1drv.ms/t/s!Ap9xAxIuzM0xlLlbUgqNsGvbe4Eq6w

#include <string>
#include <iostream>
#include <vector>
#include <fstream>
#include <sstream>
#include <map>

int main() {
    string path1 = "beam_4372-2_u.txt";
    string path2 = "beam_4372-2_d.txt";
    vector<track> vtrack1, vtrack2;

    cout << vtrack1.size() << " -> ";
    read_txt(path1, vtrack1);
    cout << vtrack1.size() << endl;

    cout << vtrack2.size() << " -> ";
    read_txt(path2, vtrack2);
    cout << vtrack2.size() << endl;

    map<string, int> hist;
    map<string, int> hist_bg;

    double bin_width = 0.002; //bin width
    int n_bin_bg = 3; //number of bin for background

    //add displacement to bins
    for (int i1 = 0; i1 < vtrack1.size(); i1++) {
        for (int i2 = 0; i2 < vtrack2.size(); i2++) {
            double dx = vtrack2[i2].px - vtrack1[i1].px;
            double dy = vtrack2[i2].py - vtrack1[i1].py;

            int dxi = round(dx / bin_width);
            int dyi = round(dy / bin_width);
            hist[to_string(dxi) + "_" + to_string(dyi)]++;
            for (int bgx = -(n_bin_bg - 1) / 2; bgx <= (n_bin_bg - 1) / 2; bgx++) {
                for (int bgy = -(n_bin_bg - 1) / 2; bgy <= (n_bin_bg - 1) / 2; bgy++) {
                    hist_bg[to_string(dxi + bgx) + "_" + to_string(dyi + bgy)]++;
                }
            }
        }
    }

    //search peak
    int max_value = 0;
    string peak_position = "";
    for (auto &p : hist) {
        p.second -= hist_bg[p.first] / double(n_bin_bg*n_bin_bg);
        if (max_value < p.second)
        {
            peak_position = p.first;
            max_value = p.second;
        }
    }

    //convert to position
    string buf;
    stringstream ss(peak_position);
    getline(ss, buf, '_');
    int peak_position_x = std::stoi(buf);
    getline(ss, buf, '_');
    int peak_posision_y = std::stoi(buf);

    //draw 7x7 bins
    for (int ix = -3; ix <= 3; ix++) {
        for (int iy = -3; iy <= 3; iy++) {
            string str = std::to_string(ix + peak_position_x) + "_" + std::to_string(iy + peak_posision_y);
            cout << hist[str] << " ";
        }
        cout << endl;
    }
    cout << "peak (x, y) = (" << peak_position_x*bin_width << ", " << peak_posision_y*bin_width << ") " << endl;
    cout << "peak bin content = " << max_value << endl;
    cin.get(); //to stop the console.
}