cctbx.xfel 読解メモ

cctbx.xfel のコードを読んでいく。メモなのでぐちゃぐちゃ。

cctbx.xfel は LCLS のシステムと強く結合している。LCLS では、生データは XTC というストリームで供給される。これを処理するのが pyana というライブラリ。XTC はイベントベースのストリームで、「設定が変わった」「XFEL パルスが来た」といったイベントがずらずらと記述されている。pyana ライブラリは、XTC ファイルを解釈し、イベントに応じて指定されたコールバック関数を呼び出す。例えば、パルスの統計値を表示する xfel/cxi/cspad_ana/mod_daq_status.py では event() でイベントのメタデータを取得して可視化している。event() は、L1Accept (level 1 filter で accept という意味か) イベントに対応するコールバック。他のコールバック関数には、beginjob, endjob, begincalibcycle などがある模様。これらをまとめて「モジュール」と呼び、xfel/cxi/cspad_ana/mod_*.py にいろいろな実装がある。モジュールは複数登録しておくことができ、hit find して、当たりを指数付けするといった流れを構築できる。ここまでで分かるように、メインループは pyana 側にあって、cctbx.xfel が従であることに注意。
# pyana と psana の関係は?

xfel/cxi/cspad_ana/mod_mar.py は XPP 実験ホールにある Rayonix detector による goniometer-based serial crystallography のためのモジュール。この場合、XTC ストリームには画像そのものは含まれず、メタデータのみ入っている模様。そこで、dxtbx を用いて、ファイルシステムから MarCCD 形式の画像を読み込み、event 構造体にセットする。以下の流れは CSPAD の場合と同じで、common_mode モジュールに引き継がれる。

xfel/cxi/cspad_ana/mod_hitfind.py はもっともコアなモジュール。xfel/cxi/cspad_ana/common_mode.py
親クラスにしているので、そちらから読む。これは、CSPAD に common mode correction などを適用する。はっきり言って、検出器ドライバの仕事である。また、cspad_tbx の image() を呼び出し、タイルを物理的な位置に再構築する。最終的な画像は、 cspad_tbx の dpack() によってビームセンタなどのメタデータと統合された Dict になる。Rayonix の場合 (init で device 引数が "marccd") は#TODO:未確認。

self.sections は検出器のパネル構成を保持している。

self.sections = calib2sections(cspad_tbx.getOptString(calib_dir))

として作られる。この関数は、xfel/cxi/cspad_ana/parse_calib.py にある。

xfel/cxi/cspad_ana/cspad_tbx.py は CSPAD のタイルを物理的な配置に再構築するルーチン image() を含む正念場。

二次元配列に再構成してしまうと、検出器パネル間の gap 領域も配列に含まれてしまう。これを除外するために、active_areas という変数が用意されている(dpack された Dict では ACTIVE_AREAS になる)。

CSPAD は 4 つの quadrant から成る。1つの quadrant は 8 つの sensor (section) からできている。1つの section は 2 つの ASIC からできている。したがって、全部で 64 のパネルがあることになる。1つの ASIC は 185x194px である。これが 4 px のギャップを挟んで 185x392px の section を作る。cbcaa() 関数に説明あり。sections の配列は、 865, 1121, 1059, 1306, 1062, 1121, 1256, 1306, ... となっているが、左上 slow, fast、右下 slow fast 4 要素 x 64 ASIC で 256 要素ある(これはパネルを再構成したあとの座標!)

これらの準備が整ったら、xfel/cxi/integrate_image_api.py が呼び出される。そして、実際の処理はxfel/cxi/display_spots.py の run_one_index() に丸投げである。

途中で xfel/detector_formats.py が参照されるが、本質的なことは何もしていない。CSPAD のジオメトリが、実験時期によって何通りかあるので、それを判定してるだけ。detector_format_version() がそれで、address (XTC イベント中での画像データの論理パス)と timestamp から判定して "CXI 10.2" のような文字列を返す。

xfel/cxi/display_spots.py は単独で実行すると、画像ビューアとして動作する。そのため display_spot という名前になっているのだ。最初の重要な仕事は、iotbx.detectors.ImageFactory を、独自のものに差し替えること。これにより、iotbx/detectors/npy.py にある NpyImage クラスを使って、quadrant と sensor の translation が適用されるようになる。指数付け・積分として動作するときの流れは、run_one_index() → run_one_index_core() である。ここからは、rstbx/new_horizons/index.py に処理が引き継がれる。ちなみに、new horizons というのは、NIH のグラント名 Realizing New Horizons in X-ray Crystallography Data Processing に由来する。。

xfel/phil_preferences.py の load_cxi_phil は、ハードコードされたデフォルト設定を読み込んだあと、設定ファイル phil をパースして適用する。ここで重要になってくるのは、params.distl.tile_translations と params.distl.quad_translations パラメータである。これが与えられていない場合、デフォルトの値として spotfinder/applications/xfel/cxi_phil.py にハードコードされた値(当然、CXI バージョンによって異なる)が利用される。後者は ASIC ごとの補正値であり、後者は 4つの quadrant ごとの補正値である。どちらも並進成分のみで回転を含まない。詳細は cctbx.xfel wiki の Metrology refinement を参照。求め方の実例は Thermolysin のチュートリアル を参照。64 の ASIC に対する、回転を含む微調節は distl ではなくて integration.subpixel_joint_model である。また、spotfinder の検出パラメータのデフォルトも定義されている。なお、驚くべきことに、

elif cxi_version in ["Sacla.MPCCD"]:

という場合分けだけは存在する。中身はないが。

rstbx/new_horizons/index.py の pre_indexing_validation と pack_names は引数を整理しているだけであまり意味はない。still でない場合(pseudo-rotation?) に意味があるのかも。new_horizons_state は、指数付けと積分に関する情報をまとめたクラスである。コンストラクタは、horizons_phil.spotfinder の値を見て、スポット検出に DISTL/LABELIT を使うか SPECK を使うか決めて、これらをインスタンス化する(が、後者も 結局 LABELIT を呼んでいるっぽい?)。process() で、実際のエンジンを作動させる。

xfel/metrology/mark10.py が最新の metrology まわりか。mark1 とか 2 とかいろいろあるが、それぞれの違いは 920 行目あたりから始まるコメントに説明がある。mosaicity といった文字も見えるが……?

xfel/command_line/cxi_view.pyxfel/cxi/display_spots.py を呼ぶだけ。実際のビューアは rstbx/viewer/ 以下にある。rstbx/viewer/frame.py でフレームを作り、画像データを rstbx/viewer/__init__.py にある image クラスとして展開する。画像データは、ファイル名としても、データそのものとしても渡せることに注意(仮引数の命名が悪くて、ファイル名を渡すように見えるけれども)。読んだ画像は image オブジェクトの _raw プロパティとして保存される。これを、コントラスト調節などして表示用の二次元画像に変換するのが create_flex_image()。実装は iotbx/detectors/display.h である。調節する UI は rstbx/viewer/calibration_frame.py にある。

#再構築されたイメージで、subpixel_jointmodel はどこでどう適用されるのか?

xfel/metrology_ext.cpp 何かやってる

#続く

以上より、新規の segmented detector に対応させるために必要な作業は、次のようになると思われる。

  • pyana に依存しない駆動系を作る。mod_hitfinder.py でやっていることを真似して データの Dict を作り、display_spots.py の run_one_index() を呼ぶ。
  • 検出器のファイル形式をサポートする。active area をきちんと定義する
  • panel metrology correction を定義する
  • metrology refinement を実装する。mark10.py を改造することになる。

はじめの2つは容易だが、残り2つは複雑。DIALS ベースにしたほうが簡単だろう。