Table of Contents... Hide

Development of MuFastSteering algorithm with NSW by T.Kumaoka

概要

開発した環境は以下のリンクのgitから持って来れます。 https://gitlab.cern.ch/tkumaoka/nswDevelop

本開発環境の扱い方は以下の画像のようになっています。

基本的にはRDOもしくはESDの生成の際に作られるntuple fileを元にし、順番に動かしその都度生成されるroot fileを次の段階で読み込ませることで最終的な性能評価がおこなえます。矢印は生成したroot fileを次に読み込ませる先を示しており、点線は点線の内のどちらかを読み込ませることを意味します。これは各ディレクトリのmain/~.hとcompile.cxxを書き換えることで決めて下さい。この開発環境は本来athenaのTriglL2MuonSAの MuFastSteeringで使われることが期待されるコードをpirvateで動かすためのものです。そのためL1や他の検出器の情報などは使えませんので注意して下さい。またほぼ原型を留めていないですが、一応基となったのは若宮さんの 修士の研究となっています。MDTとのmatching等は若宮さんが行われていたのでそちらを参考にすると良いかもしれません。

注意：

• 熊岡の修論は少し古いコードを元に書いているためあまり参考にしないでください。
• プログラム初心者であるため、用語の誤りやコードもあまり良いものでは無いかもしれませんので、もし引き継いて下さる方がいらした場合は修正して頂けると幸いです。
• 極力ここに説明は書きますが基本はコードを確認して下さい。
• 開発途中に上の方針は扱うデータ仕様が頻繁に変わったため、一度開発したものも全ては引き継げず、正直ぐちゃぐちゃな感じで終わってしまっています。あくまで性能評価の研究であって、実装するためのコードはできる人ならご自身で開発した方が早いかもしれません。

ひょっとすると返信があるかもしれないので分からないことがあればここにメールを下さい(develophltnswあっとgmail.com).

共通内容

・各ディレクトリに入り、そこで

g++ compile.cxx -o Run `root-config --cflags --libs` 

でコンパイルし、 ./Run で実行するとそれぞれで出力ファイルが生成されます。bashrcにを以下のエイリアスを作っておくとrpp compile.cxxで実行までやってくれて楽だと思います。

function rpp(){
g++ $1 -o Run `root-config --cflags --libs` && ./Run
}

・main関数はcompile.cxxで定義されていますが、基本は各ディレクトリにあるmainの下のrootの MakeClass(生成方法はリンクのp13)で生成した~.hと~.cxxがmainの枠組みです。TrigL2MuonSAであればmain/MTrigL2MuonSA.cxxがMuFastSteering.cxxの役割を担っていると考えて下さい。また読み込ませるroot file(ntuple)はmain/~.hとcompile.cxxで確認して下さい(各ディレクトリでrootと検索すると楽です)。入力fileを変更するときはこの両方を変更するようお忘れなく。また出力ファイルはmain/~.cxxで確認して下さい。基本的には各段階ごとにntuple fileを生成し、次のステップでそのfileを読み込ませることになっており、複数fileが生成されてしまうので注意して下さい。

・TrigL2MuonSAには Checkerというクラスが作ってあり、どのディレクトリでも使えるように読み込ませております。中身は TrigL2MuonSA/src/Checker.cxxを確認して下さい。基本は複数のパラメータを確認する時にcoutで書くのが面倒であるため楽にしたものです。あまりよろしいやり方ではありませんが、各main/~.hにTrigL2MuonSA::Checker c;としてグローバルなオブジェクトを作っており、どの関数内でも使えるようになっています。エラーの位置を確認したいときは、c.p0(int);とするとcheck位置がわかりやすいものとなります。また、確認したいパラメータが3つのときはc.p3("p0, p1, p2", p0, p1, p2);とすると、出力は(p0, p1, p2) = (value of p0, value of p1, value of p2) となります。p6まであるので6つまで確認できます。増やしたければChecker.hとChecker.cxxで定義すれば簡単に増やせます。(出力は全てdoubleなので注意)。実装する際にはこれは消してもらえると幸いです(c.pで検索をかければ何処にあるかはすぐに分かると思います)。

・特定のeventを確認するための方法が用意されています。Checker同様各main/~.hにグローバルな変数として、

vector<int> targetEventFilter, targetEvent;
int jentryTarget = 0;
bool targetEventTrig = 0;

が定義されており、src/の関数の好きなところで(例えばif(dEta > 0.1))targetEventTrig = 1;とするとLoop関数の最後でそのイベントの番号をtargetEventに詰めて、全てのイベントの解析が完了後event番号を吐き出してくれます。その出力をそのままevent loop(main/~.cxx)のfor文の真上あたり

//target Event Number input KumaKumaKumaKumaKumaKumaKumaKumaKumaKumaKuma
//ex targetEventFilter.push_back(#);
    targetEventFilter.push_back(0);
//target Event Number input KumaKumaKumaKumaKumaKumaKumaKumaKumaKumaKuma

にcopy pastし、

    //Event Filter KumaKumaKumaKumaKumaKumaKumaKumaKumaKumaKumaKumaKumaKuma
    targetEventTrig = 0;
    if(0){
      if(jentryTarget == targetEventFilter.size()) jentryTarget -= 1;
      if(targetEventFilter.at(jentryTarget) != jentry) continue;
      else jentryTarget++;
    }
    //Event Filter KumaKumaKumaKumaKumaKumaKumaKumaKumaKumaKumaKumaKumaKuma

のif(0)をif(1)とすると指定のイベントだけを動かしてくれるようになります。各ディレクトリ共通で使えるので、TrigL2MuonSAで出力したeventをclusteringではどのようになっているのも確かめることができます。個人的には便利だったのでよければ使ってみて下さい。

・ひとまず本コードを動かすために手打ちのNSWの値等はmain/~.hの読み込むntupleの変数の定義の下、コンストラクタの定義の上あたりに書いてあります。またそのディレクトリで使う関数等もここで定義されています。

各ディレクトリの説明

truthSegNtuple

(RDOのntuple.rootからHitsの値としてあるtruthの情報だけを引き出し扱い易いよう整形するためのpackage)

athenaとpileupのevent fileにはtruthの情報は含まれないため一度Digits_~の情報と分離する必要があるため作られたpackage。結局athenaの出力は実装担当に作って貰えていないため、athenaの性能確認としては使えていない。

clustering

(RDOのntuple.rootからDigitsの値を用いて熊岡が開発した方法でクラスタリングをし、その他の値もTrigL2MuonSAで扱い易いように整形するpackage)

コードは以前使っていた環境のものを無理やり新しい環境で使えるようにしているため、非効率なことをしていたり無駄な変数があったりします(すみません、直しきれませんでした。)。また、正直offline groupのクラスタリングを使う可能性が大きいのであまり意味のないものとなるはずです。さらにはpileup込みのサンプルにおいてはMMにバグ（layer数が16以上)あるので使えません。ただ一応ざっくりとした説明は残しておきます。

また本当は時間も測るべきでしたがそこまでできておらず、offlineのクラスタリングとどちらが良いのかは決めかねる状態です。

MmMuTpcClustering.cxxはストリップのレスポンス時間(電子のドリフト時間)の差分とストリップの位置の差分からストリップがあるzでのミューオンの入射位置を求めます(右図)。もう少し正しくいうと、もっともレスポンス時間の早かったストリップの位置に(メッシュとストリップの間隔)×(時間の差分と位置の差分から求めた傾き)から得られるdRを足します。ただし、反応すべきもっともヒット位置に近いストリップが反応しなかった時、ストリップ1本分以上の誤差となります。

StgcDataSort.cxxはpad, strip, wireを仕分け、また扱いやすいように整形する関数です。

StgcTimeMerge.cxxは同じチャンネルがdead timeよりも短い間隔で複数回データを残していることがあったため(少なくともこのコードを作った時は)、もっともレスポンス時間の短いデータだけを取ってくる関数です。

StgcChargeGausClustering.cxxは連続したストリップの電荷情報と位置情報を使ってガウスフィットし、その頂点をwireのz位置でのミューオンのRを求める関数です（左図）。

ストリップおよびクラスターの分解能とefficiencyはさらに下のディレクトリーにある effEstimateで調べることができます。

これを書いている時点ではMM(strip, 開発したもの、offline全て)のefficiencyが不自然に低いところがありますが、原因は見つかっておりません。

transePrdToCluData

(ESDのntuple.rootからoffiline groupの手法クラスタリングされた値であるPRDの値をTrigL2MuonSAで扱い易いように整形するpackage)

TrigL2MuonSA (本研究の本体)

本packageでも実装するためのクラスは処理はStgcRecoFitter.h, MmRecoFitter.h, MmStgcMerge.h、データの構造は StgcData.h, MmData.h, NswSegment.hだけである。 またLeastSquares.hは、共通して使う関数(最小二乗法の計算)であるためクラスとして定義されており、既存のpackage等を使わない場合はこれも実装する必要がある。その他の関数などは本環境にてこれらを動かすためのものである。熊岡の修士の研究では暫定的な性能評価を行うため、pTを求めることを行なっていたが、本環境ではpTは計算できない。まずこの環境に今までのpT計算のコードを組み込むより実装してしまった方が早い。また頑張って移したとしても、MDTのデータを取ってくることが難しく(以前行なっていたesdからmdtのntupleを作る方法が仕様によりそのままではできない様になっていた。）、NSWのntupleを生成するコードに書き加えてMDTを追加することもできるとは思うが、結局実装した方が早いと考えられる。

研究のはじめの方ではsTGCではpad、MMではμTPC法で計算した各層で求まる荷電粒子の入射角等を用いて再構成に使うhit情報を絞り、部分飛跡の再構成に掛かる時間の短縮するアルゴリズムを開発していた。しかし、これらのパラメータはofflineのクラスタリング後のデータ(PRD)には含まれておらず、RDOとPRD両方に対応するアルゴリズムの開発を迫られたため新しい環境には移行していない。一方で最後の最後に結局padを含む結果となったが、移行するほどの時間が無かったため使うとしても位置からの開発となると思われる。領域を絞らなくても一応時間は問題無い可能性があるため、padを用いるのphiの計算のためだけでいいかもしれない。一応指針だけは示して置こうと思う。stripで再構成したprojection（セクター中心に投影した）部分飛跡から各層に置ける通過位置を求め、そこにあるpadのchannelの情報を持ってくる。R方向を絞るだけでも8層のpadのchannelはほとんど一意に決まるはず...定めた8層のpad channelから絞られるvirtual padのphiを採用するという流れである。padの大きさがそれぞれ違うため、求めたい値がずれるという問題がある( リンク参照)。

main/MTrigL2MuonSA.cxx

1. fileのincludeが書いてある。
2. loop関数のスタート
3. 書き込むroot fileの定義(変数やbranchの設定)
4. 興味のあるeventのみをfilterするための設定
5. event loopのスタート
6. 二つ上の項目において設定したevent以外はcontinueする。
7. 想定されるdigitizationのdata構造に,eventの値を詰める( transeNtupleToHits.cxx)。~Hitは同じ情報のセット(ex. r, z, phi, etc...)をvectorにして纏めたもの。(これらはofflineの実装されているものを参考にしたが、offlineのものも中身が変わることはよくあるため定期的にチェックが必要)。
8. 1 eventで再構成するsegmentを入れていくvectorの定義:std::vector<TrigL2MuonSA::nswSegment> stgcSegmentV;
9. 本環境においてL1は存在しないため、実装の際のRoiの情報で必要となるeta, phiの情報を取ってくる関数。ただしeta, phiの値はMC truth muonのNSW検出器中心に置けるeta, phiであるため、絶対的に正しい値であり、またfake triggerなどは絶対に含まれないしようとなってしまっている。本当はeta, phiはある程度スメアしたり、fake trigger等も混ぜるべきではある。
10. make roiで作ったroiの数(MC truth muonの数)だけ回すloopのスタート.
11. 実装するための関数のスタート
12. findMmSuperPoints (実装対象) : MM検出器単体によってRoIに対し一つの部分飛跡(surper point : eta, phi, theta)を求める。
13. findStgcSuperPoints (実装対象) : sTGC検出器単体によってRoIに対し一つの部分飛跡(surper point : eta, phi, theta)を求める。
14. mmStgcMerge (実装対象) : MM検出器とsTGC検出器それぞれで再構成した部分飛跡を合わせ、NSW検出器としての単体によってRoIに対し一つの部分飛跡(surper point : eta, phi, theta)を求める。
15. event loopの終わり.
16. 出力ファイルに1eventで再構成したsegmentの情報をfillする( inputNtuple.cxx).
17. 気になったeventがあった場合、そのevent番号を記録する。if(targetEventTrig) targetEvent.push_back(jentry);
18. loop関数の終わり、root rileの書き出し.

idealTrack

(TrigL2MuonSA が理想的に動いた時、想定される最高性能を評価するためのpackage)

最後まで真面目に開発しきることができなかったため詳細は残しませんが、興味があればコードを解読もしくは直接聞いて下さい。

mcCompare

(開発したalgorithの性能評価を行うpackage)

基本的にはhistgramによって出力を確認します。

histgramの定義はあまりいい方法ではありませんが、global変数としてmain/mcCompare.hに定義されています。ここでは型だけが定義されており、グラフの詳細の定義はsrc/にある各検出器ごとのhistDef~.cxxのfileに書かれています。bin数等を変更したい場合はここで変更してください。fillする部分はまた検出器ごとに分かれており、mcCompare~.cxxに書かれています。画像出力はhistDraw~.cxxに書かれています。ある時からhistgramが多すぎると途中で落ちるという悲しい仕様になったため、writeしたら消すようにしています。またlayerごとの確認も本来できるのですが、現状offになっていると思います。確認したい場合はそれだけで動かすことを推奨します。また本当に重要なparameter(eta, phi, theta等)だけはここでは消さず、src/histWrite.cxxでroot fileに書き出される仕様になっています。また出力ファイルの名前を書き換える際はhistWrite.cxxを書き換えてください。

eventDisplayCheck

(気になったeventがどのようになっているのかを可視化するためのpackage)

これだけは普通にroot -l easyEventDisplay.cxxで動かして下さい。

用いるNtuple fileについて

はじめに読み込ませるfileは、official sampleのntupleか自分で生成したntupleを読み込ませる。RDOの段階でのntupleの値Digits_~はclusteringが行われていないため、clusteringのpackageに読み込ませる必要がある。一方ESDの段階でのntupleの値PRD_~はoffline groupのクラスタリングが行われた値である。この値を用いる場合はtransePrdToCluDataに読み困る必要がある。また、HITsは検出器に入射したtruthの情報である。これは基本RDOのntupleに含まれているが、新しいバージョンだとESDのntupleに含まれている場合がある。RDOのsTGCとMMは形式が違うため注意が必要（このroot fileを使っている他の人も困っているようでその内修正されるかもしれない。その場合はclusteringのpackageを書き直す必要がある。）。オフライングループもまだ安定していないため出力が変わることが頻繁にあり、使えない変数が出てきたり、有用な変数ができたりしているので適宜確認してもらえるとありがたい。

pileupを含んだroot fileには容量の問題からtruthの情報(Hits)が含まれていない。Hitsはtruth particleの検出器位置(各layer)での情報が並列(vectorやstructureではない)で入っています。そのため、生成時のtruth muonとは多少異なるeta phiを示します(散乱等で多少曲がるため)。最低でもmuonのtruthの情報は必要となるため、pileupを含めていないsignalのmuonのサンプルを別に作る必要がある。またいつからか理由はわかりませんが、仕様によりMMのHitsにpTを計算するために必要なhits_kineticEnergyが消滅いたしました... 本当はpTの情報がそのままあればいいのですが、あっても運動エネルギーから求めるしかありません。したがってMMでのtruthのpTの情報が欲しい時は他のパラメータが一致するsTGCのtruthから持ってくるしかありません。

    double truthRTheta = atan(truthRSlope);
    double truthKineticE = truthHits->at(hitNum).kineticE;
    double truthPt = sqrt(pow(truthKineticE * sin(truthRTheta),2) + 2*truthKineticE * truthMass) / 1000;//GeV

2020年3月の段階にてpileup sampleではMMのDigitデータにバグがあり、現状PRDを使わざるを得ない。

NSWValTree->Draw("Hits_sTGC_hitGlobalPositionY:Hits_sTGC_hitGlobalPositionX","Hits_sTGC_off_stationPhi==1&&Hits_sTGC_off_channel_type==1")

NSWValTree->Draw("Hits_MM_hitGlobalPositionY:Hits_MM_hitGlobalPositionX","Hits_MM_off_stationPhi==1")

NSWValTree->Draw("sqrt(Digits_sTGC_globalPosX*Digits_sTGC_globalPosX+Digits_sTGC_globalPosY*Digits_sTGC_globalPosY):Digits_sTGC_globalPosZ","eventNumber==0&&Digits_sTGC_stationPhi==1")

NSWValTree->Scan("4*(Digits_sTGC_multiplet -1) + (Digits_sTGC_gas_gap -1):Digits_sTGC_stationPhi?-1:Digits_sTGC_globalPosZ:sqrt(Digits_sTGC_globalPosX*Digits_sTGC_globalPosX+Digits_sTGC_globalPosY*Digits_sTGC_globalPosY)")
 

NSWValTree->Scan("Hits_sTGC_particleEncoding:Hits_sTGC_off_stationPhi:Hits_sTGC_off_gas_gap -1 + 4*(Hits_sTGC_off_multiplet -1):Hits_sTGC_hitGlobalPositionZ:sqrt(Hits_sTGC_hitGlobalPositionX*Hits_sTGC_hitGlobalPositionX+Hits_sTGC_hitGlobalPositionY*Hits_sTGC_hitGlobalPositionY)")

サンプルの作り方は基本は officialのものを参照して欲しいですが、一応熊岡が使っていた環境(root fileは無く、必要だった気がしたものだけをgitに挙げた)のリンクを貼っておきます。(https://gitlab.cern.ch/tkumaoka/nsw_mu_sample_production)

一応挙げた理由としては、pileup sampleを作る際にはなぜかiceppの環境ではできず、またlxplusではbach jobが消滅したためHTCondorでjobを投げる必要があるためです。HTCondorを使う上で必要となるファイル等はここにまとまっています。ただし、貰ったもの(川出健太郎様から)をあまり良くわからず改良して動いたからいいやの状況であるため余計な部分は多くあるかと思います。ただサンプルは一応作れます。pileup sampleは非常に重いためlxplusでは多くは作れず、作るにしても出力先はSWAN_Projectにすることをお勧めします。lxplusで行う必要があるのはneutrino, low pt, high ptのpileupを混ぜる操作だけになります。

HTCondor以外の必要なファイル等もpythonPackageというディレクトリにまとめております。

athenaへの実装について

熊岡の方でも実装するために勉強していたため、分かった範囲について,また想定していた実装について一応説明します。ただ担当が名古屋大学のグループに移ってしまったため、現状の実装状況がどのようになっているかは把握できていません。詳しい内容は名古屋の研究グループの方が私より遥かに深く理解しているはずなのでそちらに聞いて見てください。

実装について

野口さんが以前 修士の研究においてCSCをMuonSAに実装しているためCSCを参考にすると非常に実装しやすいと思われます。

MuFastSteering::findMuonSignature

1. RoI の数だけloop(for (p_roi=(muonRoIs).begin(); p_roi!=(muonRoIs).end(); ++p_roi) {)が走ります。
2. データのsetとなるm_~Hitsを初期化します。
3. barrel(if)とend cap(else)かを判断します(if ( m_recMuonRoIUtils.isBarrel(*p_roi) ))。
4. seed decondingをおそらく行う(sc = m_dataPreparator->prepareData(~))。この関数の中にsTGC, MM両方を加える必要があります。mdtやcscと異なることはTGCが無いためmuonRoadは存在しないことに注意。ここにおいては tgcを参考にした方が良いかもしれません。offline groupの clustringを参考にする場合はリンクを飛んでください。またofflineのdecodeについては sTgcRdoToPrepDataToolCore.cxxをそこで用いられる（入力）は StgcIdHelperで決まります。ここにある変数ならL2MuonSAでも扱うことができるでしょう。
5. このループ中におそらくどこでも問題ないので、m_stationFitter->findSuperPoints()やm_cscsegmentmaker->FindSuperPointCsc()と同等のNSWの関数を実装する必要があります(一応本研究はこの関数となるものを開発しました)。出力先は trackPatternsのsuperPoint[chamber(i_station)==0]であるはず(関数 TrigL2MuonSA::AlphaBetaEstimate::setMCFlagを参照)。C sideではSW,　A　SideではNSWが走るようswichできるよう設定する必要があります。またNSWはMMもsTGCも同等の検出器であるためどちらか一方でも機能します。そのため、どちらかがうまくいかなかった時片方だけでも走るようにswichできるcodeを書くことを推奨します。
6. TrigL2MuonSA::AlphaBetaEstimate::setMCFlagでLUTが使われているため、A sideにおいてはNSW用のLUTが使われるように設定する必要があります(そもそものLUTの作成もする必要有り)。

残っている研究内容(2020年3月)

• コードの実装
  • AODへのアウトプット
  • 部分飛跡の再構成アルゴリズムの実装(ATLASのコーディングルールに乗っ取った書き方へリファクタリングする必要あり)
  • decord (clustering)の見直し
  • CSCによるセグフォ問題
• ヒットレートの測定
• NSW用のLUTの作り直し
• その他(色々な直し)

参考URL

• TrigL2MuonSA : https://gitlab.cern.ch/atlas/athena/blob/master/Trigger/TrigAlgorithms/TrigL2MuonSA
• HLT tutorial： https://twiki.cern.ch/twiki/bin/viewauth/Atlas/L2MuonSA
• NSW official sample : https://twiki.cern.ch/twiki/bin/viewauth/Atlas/MuonSWTestSamples
• The method of NSW sample : https://twiki.cern.ch/twiki/bin/viewauth/Atlas/NSWdigitization
• sTGC parameter book : https://twiki.cern.ch/twiki/pub/Atlas/NSWParameterBook/sTGC.pdf
• sTGC parameter of ntuple : https://twiki.cern.ch/twiki/bin/view/Atlas/DocumentationOfAnalysisNtupleBranches
• official NSW page ： https://twiki.cern.ch/twiki/bin/view/Atlas/NewSmallWheel https://twiki.cern.ch/twiki/bin/viewauth/Atlas/NSWelectronics
• NSW TDR : https://twiki.cern.ch/twiki/pub/Atlas/NewSmallWheel/NSWTDR.pdf

Meeting Group

• Muon Software meeting https://indico.cern.ch/category/696/
• NSW Software and Performance https://indico.cern.ch/category/5279/

NSW関係者

• 偉い人 : Stefano Rosati
• pileup sample : Alexandre Laurier
• MM clusetering (micro TPC method)：Patrick Scholer

Takuya Kumaoka (熊岡 卓哉) 2020

-- TakuyaKumaoka - 2020-02-25

-- TakuyaKumaoka - 2020-02-26