はじめに

このTWikiはこれからRun2の解析を開始する学生への指南書です。標準的な解析の入門については、既に優れたインストラクションがATLASのComputingチームから提供されています(例えば SoftwareTutorailSoftwareTutorialxAODAnalysisInCMake)。また日本語のチュートリアルとしては、 粒子物理コンピューティングサマースクール(ATLASソフトウェア講習会)も説明が丁寧です。このページではそれらを適宜参照し、Run2 full dataで好きに遊べるようになることを目指し、実際的な情報を伝えていきます。

なおこのページはWYSIWYGを使って無理矢理日本語で書いていますが、そのためにTWikiのいろいろな機能が制限されています。あまり見やすくないのはそのためです。ご了承ください。

参加者が見たい信号

goodle spreadsheetを用意したので、書き込んでください。

https://docs.google.com/spreadsheets/d/1MA6snSV0IGTyv76Zysng1Y_e_SjaopxgVJ7t8ahS464/edit?usp=sharing

実践編

環境設定

login.icepp上でやることを想定します。以下を打ち込んでください。あるいは~/.bashrcに書いておきましょう。

export ATLAS_LOCAL_ROOT_BASE=/cvmfs/atlas.cern.ch/repo/ATLASLocalRootBase
alias setupATLAS='source ${ATLAS_LOCAL_ROOT_BASE}/user/atlasLocalSetup.sh'
setupATLAS

とにかくまずはデータを回してみる

以下の手順を試してみましょう。現時点では何をしているか理解しなくてokです。

mkdir YOUR_WORK_DIR
cd YOUR_WORK_DIR
mkdir source build run
cd build/
asetup AthAnalysis,21.2.68,here,slc6
mv CMakeLists.txt ../source
cmake ../source
source $TestArea/*/setup.sh
cd $TestArea/../source
acmd cmake new-skeleton MyPackage
cmake --build $TestArea
cd $TestArea/../run
athena MyPackage/MyPackageAlgJobOptions.py --filesInput /gpfs/fs7001/youhei/AnaTutorial/mc16_13TeV.361107.PowhegPythia8EvtGen_AZNLOCTEQ6L1_Zmumu.merge.AOD.e3601_s3126_r9781_r9778/AOD.11865139._001783.pool.root.1 --evtMax=10

ここでYOUR_WORK_DIRは好きなpathを設定してください。asetupで設定しているAthAnalysis 21.2.68は2018年度末現在で、最新の ASG Analysis Releaseです。詳しくは ページ下の方で。filesInputで設定しているファイルパスの先には、チュートリアル用のxAODファイルが置いてあります。(xAODとはATLASで解析に使うファイルの形式の一つです。詳しくは ページ下で。)evtMaxは回すイベント数を指定しており、今回は10イベントとしています。値を-1に設定すると、ファイル内の全イベントを回してくれます。

生成したframeworkはスケルトンと呼ばれ、イベントをただ回すだけの処理をします。続いて肉付け作業をして、イベントの中身を見ていきましょう。ただしその前に注意点として、

新しくシェルを立ち上げた時

環境を設定し直す必要があります。
bashrcを書き換えていない場合は、 setupATLASを走らせてください。そのうえで以下を実行。

cd YOUR_WORK_DIR/build
asetup
source $TestArea/*/setup.sh

コードを書き換えた時

コンパイルしましょう。

cmake --build $TestArea

build オプションでパスを指定するので、コンパイルするのはどのディレクトリからでもokです。

releaseを変えるとき

更新頻度が緩やかになったとはいえ、releaseは徐々に更新されます。新しいreleaseに変更したい時は、buildディレクトリを作り直しましょう。例えば9.9.99というreleaseに変えたい時は、以下の手順を踏みます。

cd YOUR_WORK_DIR
rm -rf build
mkdir build
cd build
asetup AthAnalysis,9.9.99,here,slc6
cmake ../source/
source $TestArea/*/setup.sh
make

イベント情報をダンプする

生成したframeworkを書き換えて、xAODに保存された情報を見ていきます。まずはイベント情報を見てみましょう。
イベント情報とはrun numberやevent numberなどの、イベントの基本的な情報で、EventInfoというオブジェクトに格納されています。

sourceディレクトリ以下のMyPackage/src/MyPackageAlg.cxxを開きましょう。execute()関数に書き加えます。

const xAOD::EventInfo* ei = 0;
CHECK( evtStore()->retrieve( ei , "EventInfo" ) ); //retrieves the event info
ATH_MSG_INFO("eventNumber = " << ei->eventNumber() ); //example printout of the event number

EventInfoクラスを使うために、ファイルの上部でインクルードしましょう。

#include "xAODEventInfo/EventInfo.h"

実際にはこれらは既にコードに書かれていて、コメントアウトされている状態なので、コメントアウトを解除するだけでokです。さてこのままではインクルードしようにも、パスがわからないのでエラーが出ます。ライブラリの依存関係をMyPackage/CMakeLists.txtに教えてあげましょう。atlas_add_libraryに以下のように追記してあげます。

atlas_add_library( MyPackageLib src/*.cxx
                   PUBLIC_HEADERS MyPackage
                   INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
                   LINK_LIBRARIES ${ROOT_LIBRARIES}
                                    AthAnalysisBaseCompsLib xAODEventInfo

cmakeして、athenaコマンドで実行してみましょう。回したイベントについて、event numberがダンプされます。

EventInfo

上の例ではevent numberを表示させただけですが、EventInfoオブジェクトは他にもいくつかのイベントレベルの情報を持っています。詳細が知りたい時は、コードを見てみましょう。Gitの EventInfoのページを見てみるか、 LXRで検索してみましょう。

ATH_MSG_INFO

上のコードでは情報のダンプはATH_MSG_INFO関数で行っています。この関数はcoutと同じように使うことができ、さらに不要な時には出力を止めさせられます。似た関数にATH_MSG_ERRORやATH_MSG_DEBUGなどがあり、jobのアウトプットレベルの設定に応じて、どのレベルのメッセージを表示させるかをコントロールできます。詳しくは このページを参照。

xAODの中身のオブジェクトを見てみる

無事イベントごとの情報が見えるようになりました。ここで一度xAODの構造を見ておきましょう。普通のrootファイルに入っているtreeのように、TFileで開いて見ることもできます。treeは何種類か用意されていますが、CollectionTreeと名付けられたtreeに主だった情報が入っています。しかしまずは予め用意されたスクリプトにより、xAODに入っているコンテナを確認してみます。

checkxAOD.py <PATH_TO_xAOD>

コンテナの名前に加えて、それぞれどの程度のサイズを食っているかが確認できます。詳細を見たいコンテナが決まったら、そのコンテナに入っている変数を見てみましょう。

root -l <PATH_TO_xAOD>
root [1] CollectionTree->Print("InDetTrackParticles*")

それぞれのコンテナについてもっと詳しい情報が知りたい時は、 下に書いてあるやり方を試してみましょう。
さてチュートリアル用のサンプルには、xAODをさらに加工したDxAODも用意してあります。DxAODは解析の種類ごとに用意されていて、それぞれの解析に合った方法でサイズを小さくしています。ここではStandard Modelのグループが用意した、2 lepton解析用のDxAODを見てみましょう。checkxAOD.pyで次のファイルを読んでみます。

/gpfs/fs7001/youhei/AnaTutorial/mc16_13TeV.361107.PowhegPythia8EvtGen_AZNLOCTEQ6L1_Zmumu.deriv.DAOD_STDM3.e3601_s3126_r10201_p3526/DAOD_STDM3.13992482._000240.pool.root.1

先ほどのxAODとの違いを確認しましょう。

HistogramやTreeを出力する

Event loopの結果をヒストグラムやtreeに出力してみましょう。EventInfoを使えるようになったので、イベントレベルの情報を吐き出しましょう。

MyPackage /src/MyPackageAlg.hでprivateメンバにヒストグラムとtreeおよびtreeのbranchに入れる変数を追加します。

TH1D* m_histAverageIntPerXing = 0;
TTree* m_tree = 0;

uint32_t m_runNumber = 0;
unsigned long long m_eventNumber = 0;
uint32_t m_lumiBlock = 0;
float m_averageIntPerXing = 0;

MyPackage /src/MyPackageAlg.cxxのinitialize()関数で生成し、histogram serviceに登録してあげます。

m_histAverageIntPerXing = new TH1D("AverageIntPerXing","AverageIntPerXing",100,0,100);
CHECK( histSvc()->regHist("/MYSTREAM/AverageIntPerXing", m_histAverageIntPerXing) );
m_tree = new TTree("tree","tree");
CHECK( histSvc()->regTree("/MYSTREAM/tree", m_tree) );
m_tree->Branch("runNumber", &m_runNumber);
m_tree->Branch("eventNumber", &m_eventNumber);
m_tree->Branch("lumiBlock", &m_lumiBlock);
m_tree->Branch("averageIntPerXing", &m_averageIntPerXing);

続いてexecute()関数内でヒストグラムとtreeに値を詰めていきます。

m_runNumber = ei->runNumber();
m_eventNumber = ei->eventNumber();
m_lumiBlock = ei->lumiBlock();
m_averageIntPerXing = ei->averageInteractionsPerCrossing();
m_histAverageIntPerXing->Fill(m_averageIntPerXing);
m_tree->Fill();

cmakeをしてathenaでコマンドで実行したらmyfile.rootというファイルができているので、rootで開いて中身を確認しましょう。

トリガーをかける

イベントを選別していきましょう。まずはトリガーを要求します。要求したいトリガーを選んだら、TrigDecisionToolと呼ばれるツールで、イベントがトリガーをパスしたか見ていきます。今回はZmumuサンプルなので、ミューオントリガーを使ってみます。

MyPackageAlg.hにTrigDecisionToolをメンバとして追加します。

asg::AnaToolHandle<Trig::TrigDecisionTool> m_tdt;

インクルードする行も追加しておきましょう。

#include "AsgTools/AnaToolHandle.h"
#include "TrigDecisionTool/TrigDecisionTool.h"

MyPackageAlg.cxxのinitialize()で、TrigDecisionToolを使う準備をします。

m_tdt.setTypeAndName("Trig::TrigDecisionTool/TrigDecisionTool");
CHECK( m_tdt.initialize() );

続いてexecute()でトリガーを判定しましょう。

bool pass_HLT_mu26_ivarmedium = m_tdt->isPassed("HLT_mu26_ivarmedium");
bool pass_HLT_mu50 = m_tdt->isPassed("HLT_mu50")

TrigDecisionTool のライブラリを忘れずにMyPackage/CMakeLists.txtに追加しておきます。

atlas_add_library( MyPackageLib src/*.cxx
                   PUBLIC_HEADERS MyPackage
                   INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
                   LINK_LIBRARIES ${ROOT_LIBRARIES}
                                    AthAnalysisBaseCompsLib xAODEventInfo TrigDecisionToolLib

athenaコマンドを実行し、イベントによってトリガーが鳴ることを確認しましょう。

トリガーの選び方

トリガーには大きく分けて、prescaleがかかっているトリガーとかかっていないトリガーがあります。前者はレートを下げるために、トリガー条件をパスしてもイベントを保存しないことがあります。一般的にlow pTやlow massのイベントを見たいのでなければ、prescaleのかかっていないトリガーを使用するほうが、統計を稼げます。どんなトリガーが自分の解析に使えるかは、 このページを見て調べてみましょう。

検出器の状態が悪いイベントを落とす

実データでは、検出器の状態が完全でないイベントが保存されていることがあります。検出器の状態はフラグで保存されているため、そのようなイベントを落としてあげましょう。

// reject event if: 
bool isMC = ei->eventType( xAOD::EventInfo::IS_SIMULATION );
bool cleanEvent = true;
if(!isMC){
  if( (ei->errorState(xAOD::EventInfo::LAr)==xAOD::EventInfo::Error ) ||
      (ei->errorState(xAOD::EventInfo::Tile)==xAOD::EventInfo::Error ) ||
      (ei->errorState(xAOD::EventInfo::SCT)==xAOD::EventInfo::Error ) ||
      (ei->isEventFlagBitSet(xAOD::EventInfo::Core, 18) ) )
    {
      cleanEvent = false;
    } // end if event flags check 
} // end if the event is data

一方全ての検出器の情報が、これらの簡易的なフラグで管理されているわけではありません。実際はデータ取得後にオフラインでデータの質を一つ一つチェックして、その時各検出器が元気に動いていたかを調べます。調べた結果はlumiblockごとにGood Run Lists (GRL) と呼ばれるファイルに記録されます。したがって解析では、このGRLを用いたイベントのフィルタリングを行います。実際にやってみましょう。

まずはMyPackage/CMakeLists.txtにGRLを扱うAsgAnalysisInterfacesのパッケージを追加します。

LINK_LIBRARIES [...] AsgAnalysisInterfaces

続いてMyPackageAlg.hにGoodRunsListSelectionToolをインクルードし、メンバとして追加します。

#include "AsgAnalysisInterfaces/IGoodRunsListSelectionTool.h"
...
asg::AnaToolHandle<IGoodRunsListSelectionTool> m_grl;

MyPackageAlg.cxxのコンストラクタで、GoodRunsListSelectionToolのコンストラクタを呼び出します。

MyPackageAlg::MyPackageAlg( const std::string& name, ISvcLocator* pSvcLocator ) :
 AthAnalysisAlgorithm( name, pSvcLocator ),
...
m_grl ("GoodRunsListSelectionTool/grl", this) {

さていよいよ使用するGRLを指定します。GRLはデータの質のチェックが進むにつれて更新されますし、自分の解析がどの検出器を使うかによっても選択肢が複数あります。 このTWikiを見て使うGRLを選びましょう。
選んだGRLをinitialize()で設定します。例えば:

std::vector<std::string> vecgrl;
vecgrl.push_back("GoodRunsLists/data15_13TeV/20170619/physics_25ns_21.0.19.xml");
vecgrl.push_back("GoodRunsLists/data16_13TeV/20180129/physics_25ns_21.0.19.xml");
vecgrl.push_back("GoodRunsLists/data17_13TeV/20180619/physics_25ns_Triggerno17e33prim.xml");
vecgrl.push_back("GoodRunsLists/data18_13TeV/20190318/physics_25ns_Triggerno17e33prim.xml");
CHECK (m_grl.setProperty("GoodRunsListVec", vecgrl));
CHECK (m_grl.setProperty("PassThrough", false));
CHECK (m_grl.initialize());

さていよいよexecute()でフィルタリングを行います。

// if data check if event passes GRL 
bool passGRL = true;
if (!isMC) {
  if (!m_grl->passRunLB(*ei)) {
    passGRL = false;
  }
} // end if the event is data

シミュレーションサンプルには検出器の状態が悪いイベントは入っていないので、実データに対してathenaコマンドを走らせて、フィルタリングの結果を確認しましょう。チュートリアル用のサンプルは以下に置いてあります。

/home/youhei/fs7001/AnaTutorial/data17_13TeV.00325713.physics_Main.merge.AOD.r10260_p3399/

ミューオンを捕まえる

MyPackage/CMakeLists.txtにミューオンのパッケージを追加します。

LINK_LIBRARIES [...] xAODMuon

EventInfo の時と同様に、MyPackageAlg.cxxでインクルードします。

#include "xAODMuon/MuonContainer.h"

execute()でミューオンを見てみましょう。

//------------
// MUONS
//------------
// get muon container of interest
const xAOD::MuonContainer* muons = 0;
CHECK (evtStore()->retrieve(muons, "Muons"));
// loop over the muons in the container
for (auto muon : *muons) {
  ATH_MSG_INFO("execute(): original muon pt/eta/phi = " << ((muon)->pt() * 0.001) << " GeV/" << (muon)->eta() << "/" << (muon)->phi());
} // end for loop over muon

これでミューオンの運動量が表示されます。

さて、通常xAODに含まれるミューオンというのは、あくまで「ミューオン候補」です。偽のミューオンも多く含まれています。そのためよりミューオンらしいものだけを、identificationしてあげる必要があります。

ミューオン選別

MyPackage/CMakeLists.txtにミューオン選別のパッケージを追加します。

LINK_LIBRARIES [...] MuonAnalysisInterfacesLib

MyPackageAlg.hに選別ツールをインクルードし、メンバとして追加します。

#include <MuonAnalysisInterfaces/IMuonSelectionTool.h>
...
asg::AnaToolHandle<CP::IMuonSelectionTool> m_muonSelection;

MyPackageAlg.cxxで初期化します。まずコンストラクタでは以下のように書き換えます。

MyPackageAlg::MyPackageAlg( const std::string& name, ISvcLocator* pSvcLocator ) :
  AthAnalysisAlgorithm( name, pSvcLocator ),
...
  m_muonSelection ("CP::MuonSelectionTool", this){

続いてinitialize()で選別の条件を決めてあげます。ここではetaの範囲を2.5に、qualityをmedium以上に設定します。qualityはTight, Medium, Loose, VeryLoose, HighPt, LowPtEfficiency が存在します。

CHECK (m_muonSelection.setProperty("MaxEta", 2.5));
CHECK (m_muonSelection.setProperty("MuQuality", (int) xAOD::Muon::Quality::Medium));
CHECK (m_muonSelection.initialize());

execute()で使ってみましょう。

for (auto muon : *muons) {
...
if ( !m_muonSelection->accept(*muon) ) continue;

athenaコマンドでいくつかのミューオンが落とされることを確認しましょう。

選別に使う変数を確認

Toolでやっていることをもうちょっと見てみます。 このページにあるように、まずはミューオンのtypeによって選別をかけています。typeを見てみましょう。

const xAOD::Muon::MuonType type = muon->muonType();
switch ( type ) {
case xAOD::Muon::Combined :
  ATH_MSG_INFO("execute(): muon is Combined");
  break;
  case xAOD::Muon::MuonStandAlone :
  ATH_MSG_INFO("execute(): muon is StandAlone");
  break;
  case xAOD::Muon::SegmentTagged :
  ATH_MSG_INFO("execute(): muon is SegmentTagged");
  break;
  case xAOD::Muon::CaloTagged :
  ATH_MSG_INFO("execute(): muon is CaloTagged");
  break;
  case xAOD::Muon::SiliconAssociatedForwardMuon :
  ATH_MSG_INFO("execute(): muon is SiliconAssociatedForwardMuon");
  break;
  default:
  ATH_MSG_INFO("execute(): unknown muon type");
}

ミューオンの運動量補正とshallow copy

今度は運動量の補正を行いましょう。これはMCでのミューオンのcalibrationを、データに合わせるものです。

先ほどのミューオン選別をスキップしてMyPackage/CMakeLists.txtにMuonAnalysisInterfacesLibを追加していない場合は、ここで追加します。
MyPackageAlg.hに校正ツールをインクルードし、メンバとして追加します。

#include "MuonAnalysisInterfaces/IMuonCalibrationAndSmearingTool.h"
...
asg::AnaToolHandle<CP::IMuonCalibrationAndSmearingTool> m_muonCalibrationAndSmearingTool;

MyPackageAlg.cxxのコンストラクタで校正ツールのコンストラクタを呼び出します。

MyPackageAlg::MyPackageAlg( const std::string& name, ISvcLocator* pSvcLocator ) :
...
m_muonCalibrationAndSmearingTool ("CP::MuonCalibrationAndSmearingTool/MuonCorrectionTool",this)
{

initialize()で初期化します。

CHECK (m_muonCalibrationAndSmearingTool.initialize());

さてこの後は少し注意が必要です。retrieveしたミューオンのコンテナにはconstが付いていたことを思い出してください。ここでやりたいことは運動量を更新することですので、constが付いているのは不都合です。しかしconstを付けるのはルールなので、外すことはできません。そこで、コンテナのコピーを作って、そちらの運動量を変更してあげます。コピーにはいくつか種類がありますが、ここでは" shallow copy"を使ってみます。

まずPATInterfacesをMyPackage/CMakeLists.txtに追加します。これはshallow copyに必須ではないですが、calibrationのエラーを感知するのに使用します。

LINK_LIBRARIES [...] PATInterfaces

MyPackageAlg.cxxに必要なツールをインクルードします。

#include "PATInterfaces/CorrectionCode.h" // to check the return correction code status of tools 
#include "xAODCore/ShallowAuxContainer.h"
#include "xAODCore/ShallowCopy.h"

execute()内のミューオンループを変更しましょう。

const xAOD::MuonContainer* muons = 0;
CHECK (evtStore()->retrieve(muons, "Muons"));
auto muons_shallowCopy = xAOD::shallowCopyContainer( *muons );
std::unique_ptr<xAOD::MuonContainer> muonsSC (muons_shallowCopy.first);
std::unique_ptr<xAOD::ShallowAuxContainer> muonsAuxSC (muons_shallowCopy.second);
// loop over the muons in the container 
for (auto muonSC : *muonsSC) {
  ATH_MSG_INFO("execute(): original muon pt/eta/phi = " << (muonSC->pt() * 0.001) << " GeV/" << muonSC->eta() << "/" << muonSC->phi());
  if ( !m_muonSelection->accept(*muonSC) ) continue;
  if(m_muonCalibrationAndSmearingTool->applyCorrection(*muonSC) != CP::CorrectionCode::Ok){
    ATH_MSG_INFO ("execute(): Problem with Muon Calibration And Smearing Tool (Error or OutOfValidityRange) ");
  }
  ATH_MSG_INFO("execute(): corrected muon pt/eta/phi = " << (muonSC->pt() * 0.001) << " GeV/" << muonSC->eta() << "/" << muonSC->phi());
} // end for loop over muons

athenaコマンドを実行し、運動量が更新されることを確認しましょう。

変数の追加

shallow copyを使うと簡単にxAODオブジェクトに変数を追加することができます。

muonSC->auxdata< TYPE >( "NEWVARIABLE" ) = XXX;

例えば

muonSC->auxdata< int >( "passMedium" ) = m_muonSelection->accept(*muonSC);

とすれば、ミューオン選別の結果を保存できます。

他のxAODオブジェクトへのリンク

ミューオンというオブジェクトは LXRでの検索結果から確認できるように、運動量に加えて、電荷、author (どのアルゴリズムで再構成されたか) などの情報を持っています。しかし例えばミューオンがinner detectorでどんな飛跡を描いたのか、という情報は入っていません。これはTrackParticleというオブジェクトに含まれています。ミューオンは対応するTrackParticleへのリンクを持っているため、このリンクを辿ってTrackParticleオブジェクトを入手し、そこから必要な情報を引き出せる、という仕組みになっています。
実際にやってみましょう。ミューオンループの中に以下の行を追加します。

ElementLink< xAOD::TrackParticleContainer > trkLink;
if (muonSC->isAvailable< ElementLink< xAOD::TrackParticleContainer > > ("inDetTrackParticleLink") ){
  trkLink = muonSC->auxdata< ElementLink< xAOD::TrackParticleContainer> >("inDetTrackParticleLink");
  if(trkLink.isValid()){
    ATH_MSG_INFO("execute(): ID track pt/eta/phi/d0 = " << (*trkLink)->pt() << "/" << (*trkLink)->eta() << "/" << (*trkLink)->phi() << "/" << (*trkLink)->d0());
  }
}

athenaコマンドを実行して、trackの情報がダンプされるのを確認しましょう。注意点として、ミューオンはミューオンスペクトロメータだけで再構成されることもあります。このようなミューオンはinner trackを持ちません。この場合、リンクからはNULLが返ってきます。

ところでxAOD::Muonには trackParticle()という関数が用意されていて、リンクに直接触らずとも、TrackParticleを取ってくることもできます。これを使うと上の操作がより簡略化されます。

const xAOD::TrackParticle *IDtrack = muonSC->trackParticle(xAOD::Muon::TrackParticleType::InnerDetectorTrackParticle);
if(IDtrack) {
  ATH_MSG_INFO("execute(): ID track pt/eta/phi/d0 = " << IDtrack->pt() << "/" << IDtrack->eta() << "/" << IDtrack->phi() << "/" << IDtrack->d0());
}

やっていることは同じなので、どちらを使っても結果は変わりません。

ミューオンisolation

ミューオン選別によって本物のミューオンを取って来れはしますが、heavy flavor quarkのsemileptonic decayで出てくるnon-promptミューオンは排除できません。これを落とすためにisolationを要求してみましょう。

MyPackage /CMakeLists.txtにisolationツールのパッケージを追加します。

LINK_LIBRARIES [...] IsolationSelectionLib

MyPackageAlg.hにisolationツールをインクルードし、メンバとして追加します。

#include "IsolationSelection/IIsolationSelectionTool.h"
...
asg::AnaToolHandle<CP::IIsolationSelectionTool> m_isoTool;

MyPackageAlg.cxxで初期化します。まずコンストラクタでは以下のように書き換えます。

MyPackageAlg::MyPackageAlg( const std::string& name, ISvcLocator* pSvcLocator ) :
AthAnalysisAlgorithm( name, pSvcLocator ),
...
m_isoTool("CP::IsolationSelectionTool/IsolationSelectionTool"){

続いてinitialize()で要求するisolationの厳しさ (Working Point, WP) を決めてあげます。ここではFixedCutTightTrackOnlyを使用します。 他のWPについては isolatonのTWikiを参考にしてください。

CHECK( m_isoTool.setProperty("MuonWP", "FixedCutTightTrackOnly") );
CHECK( m_isoTool.retrieve() );

execute()で実際に使ってみましょう。

for (auto muonSC : *muonsSC) {
...
bool pass_iso = m_isoTool->accept(*muonSC);

athenaコマンドでいくつかのミューオンが落とされることを確認しましょう。

複数のミューオンを組み合わせる

例えばZボソンが2つのミューオンに崩壊した事象を見たい場合、2ミューオンで不変質量を組みたくなります。そのためにはいくつか方法がありますが、一つの素直なやり方としては、ミューオンループの結果をvectorに詰めておいて、それを後で解析することでしょう。そのためには

int m_nMuons;
std::vector<double> *m_muon_phi;
std::vector<double> *m_muon_eta;
std::vector<double> *m_muon_pt;
std::vector<double> *m_muon_e;
std::vector<double> *m_muon_charge;
std::vector<int> *m_muon_passMedium;

のような変数を追加しておけば、ループの後で改めてfor文を回して質量を計算できます。framework内でvectorのループを回してもよいですし、一度vectorをtreeに出力しておき、後から簡単なマクロで解析してもよいでしょう。
しかしここでは別のやり方を試してみましょう。framework内で質量を組みます。shallow copyを使えば運動量の更新や変数の追加が行えるので、もう一度iteratorでループを回せばよいです。

xAOD::MuonContainer::const_iterator mu_itr1 = muonsSC->begin();
xAOD::MuonContainer::const_iterator mu_end = muonsSC->end();
for ( ; mu_itr1!=mu_end; ++mu_itr1 ) {
  if(!(*mu_itr1)->auxdata< int >( "passMedium" )) continue;
  // second muon loop 
  xAOD::MuonContainer::const_iterator mu_itr2 = mu_itr1;
  for ( mu_itr2++; mu_itr2!=mu_end; ++mu_itr2 ) {
    if(!(*mu_itr2)->auxdata< int >( "passMedium" )) continue;
    if( (*mu_itr1)->charge()*(*mu_itr2)->charge() > 0 ) continue;
    // selection like if((*mu_itr1)->charge() * (*mu_itr2)->charge() > 0) continue;
    // then calculate di-muon mass
    const double dimuMass = ((*mu_itr1)->p4() + (*mu_itr2)->p4()).M();
...

ここでp4()はTLorentzVectorを返す関数で、xAOD::Muonの他、いろいろな物理オブジェクトに備えられています。

MC truth

シミュレーションではどの粒子がどんな運動量を持っているか、どこで崩壊したか、全て「答え」が分かっています。現実世界でどう見えるかをシミュレートするのではなく、直接その答えを教えてくれるのがMC truthです。

MC truthの構造

xAODのオブジェクトの中でも少し特殊なので、簡単に説明しておきます。

ATLASにおいてMC truthとは、event generatorの出力と、Geant4シミュレーションの出力を指します。ここで前者はダイアグラムレベルの情報、parton shower, hadronization, 検出器までの終状態全てを含みます。つまり検出器によるdegitization, reconstructionを除く部分は、全てMC truthのカバー範囲といえます。
xAODに含まれるMC truthは、3つのクラスで記述されます。

  • xAOD::TruthEvent
  • xAOD::TruthParticle
  • xAOD::TruthVertex
それぞれの詳細な説明は このページを見てもらうとして、ここではtruth particleとtruth vertexの関係を覚えてください。下の図のようにtruth paticleの衝突点、または崩壊点はtruth vertexで記述されるので、truth paticle -> truth vertex -> truth particle...というように、チェーンの構造になっています。

したがって(全てのチェーン構造がファイルに保存されているかはevent generator, ファイル形式にもよりますが)TruthParticleとTruthVertexを交互に辿っていくことで、陽子陽子衝突からdegitizationまでの、計算された過程を全て再現することが可能です。

さて続いて実際にMC truthを見ていきますが、2つのやり方を試してみましょう。

コンテナをretrieveしてdecay chainを辿る

parent particle -> decay vertex -> child particle -> decay vertex -> grandchild particle ...と上流からどんどん下流に辿っていきます。ここでは再帰的な関数を使ってみましょう。
まずMyPackage/CMakeLists.txtにMC truthパッケージを追加します。

LINK_LIBRARIES [...] xAODTruth

MyPackageAlg.hでTruth Particleパッケージをインクルードし、新しい関数を定義します。

#include "xAODTruth/TruthParticleContainer.h"
...
StatusCode DumpDecayChain(const xAOD::TruthParticle* parent);

MyPackageAlg.cxxでTruth Vertexパッケージをインクルードします。

#include "xAODTruth/TruthVertexContainer.h"

execute()でMC truth particleコンテナをretrieveし、その中からZボソンを探します。見つけたら先ほどの再帰的な関数を呼び出しましょう。

//------------ 
// MC Truth 
//------------ 
const xAOD::TruthParticleContainer* mctruths = 0;
CHECK ( evtStore()->retrieve( mctruths, "TruthParticles" ) );
for (auto mctruth : *mctruths) {
  if(!mctruth) continue;
  float mctruth_pt = mctruth->pt();
  if(mctruth_pt<0.00001) continue;
  if(mctruth->pdgId() == 23 && mctruth->hasDecayVtx()) {
    ATH_MSG_INFO("execute(): parent MC truth barcode" << mctruth->barcode() << " pdgId/status/pt = " << mctruth->pdgId() << "/" << mctruth->status() << "/" << mctruth_pt);
    DumpDecayChain(mctruth);
  }
}

ファイルの下部で、関数の中身を書きましょう。

StatusCode MyPackageAlg::DumpDecayChain(const xAOD::TruthParticle* parent) {
  bool hasDecayVtx = parent->hasDecayVtx();
  if(!hasDecayVtx) return StatusCode::SUCCESS;
  const xAOD::TruthVertex* vtx = parent->decayVtx();
  float vtx_x = vtx->x();
  float vtx_y = vtx->y();
  float vtx_z = vtx->z();
  float vtx_r = sqrt(vtx_x*vtx_x + vtx_y*vtx_y);
  float vtx_phi = atan2(vtx_x,vtx_y);
  ATH_MSG_INFO("DumpDecayChain(): barcode" << parent->barcode() << " decays at Z/R/phi = " << vtx_z << "/" << vtx_r << "/" << vtx_phi);
  for(unsigned int i=0; i<vtx->nOutgoingParticles(); i++) {
    const xAOD::TruthParticle* child = vtx->outgoingParticle(i);
    if(!child) continue;
    float child_pt = child->pt();
    if(child_pt<0.00001) continue;
    ATH_MSG_INFO("DumpDecayChain(): child of " << parent->barcode() << ", barcode" << child->barcode() << " pdgId/status/pt = " << child->pdgId() << "/" << child->status() << "/" << child_pt);
    DumpDecayChain(child);
  }
return StatusCode::SUCCESS;
}

athenaコマンドを実行し、Zボソンがシミュレーション中で崩壊する過程を見てみましょう。

再構成されたオブジェクトにリンクしたMC truthを取得する

MC truthのコンテナを回さなくても、再構成されたオブジェクトから直接truthのリンクを持って来ることもできます。ミューオンループ内で以下のように書いてみましょう。

bool foundlink = false;
ElementLink< xAOD::TruthParticleContainer > truthLink;
if (muonSC->isAvailable< ElementLink< xAOD::TruthParticleContainer > > ("truthParticleLink") ){
  truthLink = muonSC->auxdata< ElementLink< xAOD::TruthParticleContainer> >("truthParticleLink");
  if(truthLink.isValid()){
    foundlink = true;
  }
  if (foundlink){
    ATH_MSG_INFO("execute(): TruthParticle linked to reco muon barcode" << (*truthLink)->barcode() << " pdgId/status/pt = " << (*truthLink)->pdgId() << "/" << (*truthLink)->status() << "/" << (*truthLink)->pt());
  }
}

単純にtruthと再構成されたオブジェクトの比較がしたい場合は、こちらが便利です。

ジェットを捕まえる

まずMyPackage/CMakeLists.txtにxAODJetパッケージを追加します。

LINK_LIBRARIES [...] xAODJet

MyPackageAlg.cxxでJetコンテナパッケージをインクルードします。

#include "xAODJet/JetContainer.h"

execute()でジェットのループを回します。

// get jet container of interest 
const xAOD::JetContainer* jets = 0;
CHECK(evtStore()->retrieve( jets, "AntiKt4EMTopoJets" ));
ATH_MSG_INFO ("execute(): number of jets = " << jets->size());
// loop over the jets in the container 
for (auto jet : *jets) {
  ATH_MSG_INFO("execute(): original jet pt/eta/phi = " << (jet->pt() * 0.001) << " GeV/" << jet->eta() << "/" << jet->phi());
}

使用できるジェットの種類については ページ下で説明します。

Jet energy補正

ミューオンの場合と同様、エネルギーの補正を行います。ミューオンの場合と違うのは、実データに対しても補正を行う点です。まずはライブラリをMyPackage/CMakeLists.txtに追加します。

LINK_LIBRARIES [...] JetCalibToolsLib

続いてMyPackageAlg.hにてツールを宣言します。

#include "JetCalibTools/IJetCalibrationTool.h"
...
asg::AnaToolHandle<IJetCalibrationTool> m_jetCalibration;

MyPackageAlg.cxxのinitialize()で、補正の設定を行います。

 TString jetAlgo = "AntiKt4EMTopo"; // Jet collection, for example AntiKt4EMTopo or AntiKt4LCTopo (see below) 
 TString config = "JES_MC16Recommendation_Consolidated_EMTopo_Apr2019_Rel21.config"; // Global config (see below) 
 TString calibSeq = "JetArea_Residual_EtaJES_GSC_Smear"; // Calibration sequence to apply (see below) 
 TString calibArea = "00-04-82"; // Calibration Area tag (see below) 
 bool isData = false; // bool describing if the events are data or from simulation 
 m_jetCalibration.setTypeAndName("JetCalibrationTool/JetCalibrationTool_EMTopo_MC");
 if( !m_jetCalibration.isUserConfigured() ){
   CHECK( m_jetCalibration.setProperty("JetCollection",jetAlgo.Data()) );
   CHECK( m_jetCalibration.setProperty("ConfigFile",config.Data()) );
   CHECK( m_jetCalibration.setProperty("CalibSequence",calibSeq.Data()) );
   CHECK( m_jetCalibration.setProperty("CalibArea",calibArea.Data()) );
   CHECK( m_jetCalibration.setProperty("IsData",isData) );
   CHECK( m_jetCalibration.retrieve() );
 }

ここでjetAlgoはJetの種類を指定しています(ジェットの種類については ページ下参照)。config, calibAreaはjetAlgoに応じて変えます。
で、ここから注意が必要なのですが、calibSeqとisDataは、dataとMCに応じて変えないといけません。自動でdataかMCかを判別する方法もありますが、まずは手で書きかえる素朴なやり方を試してください。ちなみに上の例ではMCを使用します。dataの場合は:

TString calibSeq = "JetArea_Residual_EtaJES_GSC_Insitu"
bool isData = true; // bool describing if the events are data or from simulation

となります。

さらに、MCであっても、AFIIと呼ばれる高速のシミュレーションを使ったサンプルを使用する場合は、configを変えないといけません。

TString config = "JES_MC16Recommendation_AFII_EMTopo_Apr2019_Rel21.config"

ここらへんは詳しくは JetEtMissグループのTWikiを参照してください。

続いてexecute()内のジェットループを変更しましょう。

// get jet container of interest 
const xAOD::JetContainer* jets = 0;
CHECK(evtStore()->retrieve( jets, "AntiKt4EMTopoJets" ));
ATH_MSG_INFO ("execute(): number of jets = " << jets->size());
auto jets_shallowCopy = xAOD::shallowCopyContainer( *jets );
std::unique_ptr<xAOD::JetContainer> jetsSC (jets_shallowCopy.first);
std::unique_ptr<xAOD::ShallowAuxContainer> jetsAuxSC (jets_shallowCopy.second);
// loop over the jets in the container 
for (auto jetSC : *jetsSC) {
  // ATH_MSG_INFO ("execute(): jet pt = " << (jet->pt() * 0.001) << " GeV"); // just to print out something 
  ATH_MSG_INFO("execute(): original jet pt/eta/phi = " << (jetSC->pt() * 0.001) << " GeV/" << jetSC->eta() << "/" << jetSC->phi());
  m_jetCalibration->applyCalibration(*jetSC);
  ATH_MSG_INFO("execute(): corrected jet pt/eta/phi = " << (jetSC->pt() * 0.001) << " GeV/" << jetSC->eta() << "/" << jetSC->phi());
} // end for loop over jets

athenaコマンドを実行すれば、補正前後のJetの運動量が確認できます。

Jet cleaning

カロリメータのノイズや、ビームハローによって誤って再構成されたジェットを落とします。ライブラリをMyPackage/CMakeLists.txtに追加しましょう。

LINK_LIBRARIES [...] JetSelectorToolsLib

MyPackageAlg.hでインクルードします。

#include "JetSelectorTools/JetCleaningTool.h"
...
asg::AnaToolHandle<IJetSelector> m_jetCleaningTool;

MyPackageAlg.cxxにおいて、コンストラクタに以下の行を追加します(, を適宜つけること)。

MyPackageAlg::MyPackageAlg( const std::string& name, ISvcLocator* pSvcLocator ) :
AthAnalysisAlgorithm( name, pSvcLocator ),
...
m_jetCleaningTool("JetCleaningTool/JetCleaningTool")

initialize()でツールの初期化を行います。

CHECK( m_jetCleaningTool.setProperty("CutLevel", "LooseBad") );
CHECK( m_jetCleaningTool.setProperty("DoUgly", false) );
CHECK( m_jetCleaningTool.retrieve() );

execute()のジェットループにおいて、ジェットの選別を行います。例としてpass_Cleaningという変数に、選別結果を入れておきましょう。

for (auto jetSC : *jetsSC) {
...
bool pass_Cleaning = m_jetCleaningTool->keep(*jetSC);
...
}

bad jetがいるイベントでは、該当するジェット以外にも悪い影響が見られることがわかっています。そのため通常そのジェットを排除するだけでなく、イベント自体をvetoしてしまいます。詳しくは JetEtMissのTWikiを確認してください。

Pileup removal

JVTとfJVTを用いて、hard scattering由来のジェットと、pileup vertex由来のジェットを切り分けます。なおこの作業は通常METの計算前にやっておくことが求められます。

ライブラリをMyPackage/CMakeLists.txtに追加しましょう。

LINK_LIBRARIES [...] JetMomentToolsLib

MyPackageAlg.hでインクルードします。

#include "JetMomentTools/JetVertexTaggerTool.h"
#include "JetMomentTools/JetForwardJvtTool.h"
...
asg::AnaToolHandle<IJetUpdateJvt> m_jetJvtUpdateTool;
asg::AnaToolHandle<IJetModifier> m_jetFwdJvtTool;

MyPackageAlg.cxxにおいて、コンストラクタに以下の行を追加します(, を適宜つけること)。

m_jetJvtUpdateTool("JetVertexTaggerTool/JetVertexTaggerTool"),
m_jetFwdJvtTool("JetForwardJvtTool/FJVTTool")

initialize()で初期化しましょう。

CHECK( m_jetJvtUpdateTool.retrieve() );
CHECK( m_jetFwdJvtTool.setProperty("UseTightOP", true) );
CHECK( m_jetFwdJvtTool.setProperty("CentralMaxPt",60e3) );
CHECK( m_jetFwdJvtTool.setProperty("OutputDec", "passFJvt") );
CHECK( m_jetFwdJvtTool.retrieve() );

execute()で実際にツールを使用します。

for (auto jetSC : *jetsSC) {
  ...
  float jvt = m_jetJvtUpdateTool->updateJvt(*jetSC);
  ...
}
m_jetFwdJvtTool->modify(*jetsSC);

この結果、 jetのauxdecor<float>("Jvt") という変数に入れられたJVTの値が更新され、passFJvtという変数に、fJVTの結果が詰められます。

Jet Tile Collection

不要?

b-tag

Jetの種類を分類します。tauやquark-gluon taggingに加え、最近はboson taggingなんかもありますが、ここでは最も一般的なb-tagを試しましょう。まずはライブラリをMyPackage/CMakeLists.txtに追加します。

LINK_LIBRARIES [...] xAODBTaggingEfficiency

続いてMyPackageAlg.hでツールを宣言します。

#include "xAODBTaggingEfficiency/BTaggingSelectionTool.h"
...
asg::AnaToolHandle<IBTaggingSelectionTool> m_btagSelectionTool;

MyPackageAlg.cxxのコンストラクタで、BTaggingSelectionToolのコンストラクタを呼び出します。

MyPackageAlg::MyPackageAlg( const std::string& name, ISvcLocator* pSvcLocator ) :
  AthAnalysisAlgorithm( name, pSvcLocator ),
  ...
  m_btagSelectionTool ("BTaggingSelectionTool", this)

MyPackageAlg.cxxのinitialize()で初期化します。

CHECK( m_btagSelectionTool.setProperty( "FlvTagCutDefinitionsFileName","xAODBTaggingEfficiency/13TeV/2017-21-13TeV-MC16-CDI-2018-02-09_v1.root" ) );
CHECK( m_btagSelectionTool.setProperty("TaggerName", "MV2c10" ) );
CHECK( m_btagSelectionTool.setProperty("OperatingPoint", "FixedCutBEff_70") );
CHECK( m_btagSelectionTool.setProperty("JetAuthor", "AntiKt4EMTopoJets" ) );
CHECK( m_btagSelectionTool.retrieve() );

最初の行は閾値の情報が入っているCDIファイルを指定しています。閾値はb-tagの研究が進むにつれて変更されることがあるので、xAODではなく外部ファイルに格納されています。2行目はb-tagの名前を指定していて、ここではMV2アルゴリズムを用いたMV2c10 (c-jetをlight jetの中に10%混ぜたサンプルでtrainingされたMV2 score) というものを指定しています。3行目は複数存在するworking pointの一つを指定しています。ここではpTに対して一定のefficiency「ではない」けれど、ttbar事象に対してinclusiveに70 %となるworking pointを使用しています。他にどんなオプションがあるかについては、 FlavorTaggingグループのTWikiを見てください。

続いてexecute()内のジェットループを変更して、b-tagを行いましょう。

for (auto jetSC : *jetsSC) {
  ...
  bool tagged = m_btagSelectionTool->accept(*jetSC);
}

これだけでokです。b-tag scoreが欲しい場合、さらに閾値の値が欲しい場合は次のようにします。

  double tagweight = 0;
  m_btagSelectionTool->getTaggerWeight( *jetSC ,tagweight)
  double cutvalue = 0;
  m_btagSelectionTool->getCutValue(jetSC->pt(), cutvalue);

さて次にMC truthの情報を使って、本来は何のflavorなのかを見てみましょう。truthの情報を使うので、当然dataに対しては使えませんので念のため。

int truthLabel = -1;
 jetSC->getAttribute("HadronConeExclExtendedTruthLabelID",truthLabel);

この値が5, 54, 55の場合はtruthにb-jetということになります。
他の値の意味については こちらのTWikiを見てください。

Photonを捕まえる

このセクションではphotonを扱うので、使用するサンプルもそれに合わせて以下を想定します。

mc16a: mc16_13TeV.343981.PowhegPythia8EvtGen_NNLOPS_nnlo_30_ggH125_gamgam.deriv.DAOD_HIGG1D1.e5607_s3126_r9364_r9315_p3832
mc16d: mc16_13TeV.343981.PowhegPythia8EvtGen_NNLOPS_nnlo_30_ggH125_gamgam.deriv.DAOD_HIGG1D1.e5607_e5984_s3126_r10201_r10210_p3832
mc16e: mc16_13TeV.343981.PowhegPythia8EvtGen_NNLOPS_nnlo_30_ggH125_gamgam.deriv.DAOD_HIGG1D1.e5607_e5984_s3126_r10724_r10726_p3832
data15: data15_13TeV.periodAllYear.physics_Main.PhysCont.DAOD_HIGG1D1.grp15_v01_p3704
data16: data16_13TeV.periodAllYear.physics_Main.PhysCont.DAOD_HIGG1D1.grp16_v01_p3704
data17: data17_13TeV.periodAllYear.physics_Main.PhysCont.DAOD_HIGG1D1.grp17_v01_p3704
data18: data18_13TeV.periodAllYear.physics_Main.PhysCont.DAOD_HIGG1D1.grp18_v01_p3704

photonはEGammaグループの管轄なので、xAODEgammaをCMakeLists.txtに追加します。

LINK_LIBRARIES [...] xAODEgamma

続いてMyPackageAlg.cxxでコンテナパッケージをインクルードします。

#include "xAODEgamma/PhotonContainer.h"

さらにexecute()でphotonのループを回しましょう。以下では今後のために最初からshallow copyを導入しています。

//------------ 
// PHOTONS 
//------------ 
// get photon container of interest 
const xAOD::PhotonContainer* photons = 0;
CHECK (evtStore()->retrieve(photons, "Photons"));
auto photons_shallowCopy = xAOD::shallowCopyContainer( *photons );
std::unique_ptr<xAOD::PhotonContainer> photonsSC (photons_shallowCopy.first);
std::unique_ptr<xAOD::ShallowAuxContainer> photonsAuxSC (photons_shallowCopy.second);
// loop over the photons in the container 
for (auto photonSC : *photonsSC) {
  ATH_MSG_INFO("execute(): original photon pt/eta/phi = " << (photonSC->pt() * 0.001) << " GeV/" << photonSC->eta() << "/" << photonSC->phi());
} // end for loop over photons

エネルギー補正

photonの場合はdataとMCの両方に補正をかけます。まずはCMakeLists.txtにライブラリを書き加えましょう。

LINK_LIBRARIES [...] EgammaAnalysisInterfacesLib

続いてMyPackageAlg.hでツールを宣言します。

#include "EgammaAnalysisInterfaces/IEgammaCalibrationAndSmearingTool.h"
...
asg::AnaToolHandle<CP::IEgammaCalibrationAndSmearingTool> m_EgammaCalibrationAndSmearingTool;

さてMyPackageAlg.cxxで宣言したツールを使用します。まずはコンストラクタを以下のように書き換えます。

MyPackageAlg::MyPackageAlg( const std::string& name, ISvcLocator* pSvcLocator ) :
  AthAnalysisAlgorithm( name, pSvcLocator ),
  ...
  m_EgammaCalibrationAndSmearingTool("CP::EgammaCalibrationAndSmearingTool/EgammaCalibrationAndSmearingTool", this)

initialize()では、補正に使うモデルを指定してあげます。

CHECK( m_EgammaCalibrationAndSmearingTool.setProperty("ESModel", "es2018_R21_v0") );
CHECK( m_EgammaCalibrationAndSmearingTool.retrieve() );

注意1:このページで 最初に設定したrelease 21.2.68は既に少し古くなっており、2019年6月現在最新の"es2018_R21_v0"は使えません。そのままではジョブを走らせた時に怒られます。解決策の一つは ページ上部に書いてあるreleaseの変更方法を使って、新しいreleaseに更新することです。AthAnalysis,21.2.74以降なら動きます。releaseを変えない場合は、"es2018_R21_v0"を"es2017_R21_v1"に替えてください。

さてexecute()のphotonのループ内で補正を行います。

for (auto photonSC : *photonsSC) {
...
  if(m_EgammaCalibrationAndSmearingTool->applyCorrection(*photonSC) != CP::CorrectionCode::Ok){
    ATH_MSG_INFO ("execute(): Problem with Photon Calibration And Smearing Tool (Error or OutOfValidityRange)");
  }

athenaコマンドでエネルギーの補正を確認してみましょう。

注意2: MCに対して走らせた場合、pileup toolどうこうというエラーメッセージが出ます。これはMCでpileup環境がdataをちゃんと反映していないことへの警告文です。とりあえずそのまま放置しても動きますが、気になる場合は 下のpileup reweightingのインストラクションを先に読んでください。

シャワー形状補正と選別

photonの選別にはシャワー形状を使用しますが、MCとdataで形状にずれがあるので、MCに対して補正を与えます。

注意: AFIIサンプルには使用しません。

MyPackageAlg.hでツールを宣言します。

#include "EgammaAnalysisInterfaces/IElectronPhotonShowerShapeFudgeTool.h"
...
asg::AnaToolHandle<IElectronPhotonShowerShapeFudgeTool> m_fudgeMCTool;

MyPackageAlg.cxxで宣言したツールを使用します。いつもの手順でコンストラクタを書き換えます。

m_fudgeMCTool("ElectronPhotonShowerShapeFudgeTool/ElectronPhotonShowerShapeFudgeTool",this)

続いてinitialize()

CHECK( m_fudgeMCTool.setProperty("Preselection", 22) );
CHECK( m_fudgeMCTool.setProperty("FFCalibFile","ElectronPhotonShowerShapeFudgeTool/v2/PhotonFudgeFactors.root") );
CHECK( m_fudgeMCTool.retrieve() );

そしてexecute()のphotonループで補正を与えます。

if(isMC) { // NOT apply Fudge to AFII samples!!!!!!!!!!!!!!!!!!!!!!!!! 
  float unfudged_weta1 = photonSC->auxdata<float>("weta1");
  if(m_fudgeMCTool->applyCorrection(*photonSC) != CP::CorrectionCode::Ok){
    ATH_MSG_INFO ("execute(): Problem with Photon Shower Shift Correction (Error or OutOfValidityRange)");
  }
  float fudged_weta1 = photonSC->auxdata<float>("weta1");
  ATH_MSG_INFO ("execute(): weta1 = " << unfudged_weta1 << " -> " << fudged_weta1);
}

続いてphotonの選別を行っていきます。
まずはphotonのenergy depositが検出器の有効範囲に入っていることを要求します。energyは主にEMカロリメータの第2層で測定するため、この層でのeta座標を用います。

bool pass_eta = false;
if( fabs( photonSC->caloCluster()->etaBE(2) ) < 1.37
    || (fabs( photonSC->caloCluster()->etaBE(2) ) > 1.52
        && fabs( photonSC->caloCluster()->etaBE(2) ) < 2.37) ) {
  pass_eta = true;
}

ここで1.37 < eta(s2) < 1.52と、ある範囲を除いています。これはcrack regionと呼ばれる領域で、ここでfake photonのrejection等が悪化するため、今回の解析からは除きます。

次に(特殊な信号を仮定し ない限り) 明らかに異常なphotonを落とす手続きを行います。

bool pass_badcls = photonSC->isGoodOQ(xAOD::EgammaParameters::BADCLUSPHOTON);
bool pass_OQ = false;
if(!( (photonSC->OQ()&1073741824)!=0 ||
      ( (photonSC->OQ()&134217728)!=0 &&
        (photonSC->showerShapeValue(xAOD::EgammaParameters::Reta) >0.98
         ||photonSC->showerShapeValue(xAOD::EgammaParameters::f1) > 0.4
         ||(photonSC->OQ()&67108864)!=0)
        ) ) ){
  pass_OQ = true;
}

author cutと呼ばれる手続きを行います。photonコンテナには時々jetのようにtopo clusterを使用して再構成されたphotonが入っており、ある種の解析ではこれを使用するのですが、通常は用いません。

bool pass_author = false;
uint16_t author = photonSC->author();
if ((author & xAOD::EgammaParameters::AuthorPhoton) || (author & xAOD::EgammaParameters::AuthorAmbiguous)){
  pass_author = true;
}

では本番の選別を行っていきます。ついでに2016年にあったEMカロリメータのHV問題によって影響を受けたphotonの排除も行います。
MyPackageAlg.hでツールを宣言します。

#include "EgammaAnalysisInterfaces/IAsgPhotonIsEMSelector.h"
#include "EgammaAnalysisInterfaces/IAsgDeadHVCellRemovalTool.h"
...
asg::AnaToolHandle<IAsgPhotonIsEMSelector> m_photonLooseIsEMSelector;
asg::AnaToolHandle<IAsgPhotonIsEMSelector> m_photonTightIsEMSelector;
asg::AnaToolHandle<IAsgDeadHVCellRemovalTool> m_deadHVTool;

MyPackageAlg.cxxのコンストラクタをに以下を追加します。(,を適宜つけるのを忘れないこと)

m_photonLooseIsEMSelector("AsgPhotonIsEMSelector/PhotonSelIsEM_Loose",this)
m_photonTightIsEMSelector("AsgPhotonIsEMSelector/PhotonSelIsEM_Tight",this)
m_deadHVTool("AsgDeadHVCellRemovalTool/deadHVTool",this)

initialize()ではそれぞれのqualityを設定します。

CHECK( m_photonLooseIsEMSelector.setProperty("WorkingPoint", "LoosePhoton") );
CHECK( m_photonLooseIsEMSelector.retrieve() );
CHECK( m_photonTightIsEMSelector.setProperty("WorkingPoint", "TightPhoton") );
CHECK( m_photonTightIsEMSelector.retrieve() );
CHECK( m_deadHVTool.retrieve() );

execute()でツールを走らせます。
注意: シャワー形状の補正後に行ってください。

bool pass_LooseID = m_photonLooseIsEMSelector->accept(photonSC);
bool pass_TightID = m_photonTightIsEMSelector->accept(photonSC);
bool pass_HV = m_deadHVTool->accept(photonSC);

di-photon解析におけるprimary vertex選別

以下でdi-photonチャンネルでHiggsを捕まえるのに必要な技術的な情報を載せておきます。具体的にはprimary vertexの選別と、それに伴う運動量、isolationの補正を記述します。primary vertexの選別は他のオブジェクトにも影響しますが (例えばJetのJVTやleptonのimpact prameter) 、ここでは割愛します。
di-photon解析で使用するtriggerやイベント選別については、HGamグループのノートを参考にしてください。
まずはCMakeLists.txtの書き換えです。

LINK_LIBRARIES [...] PhotonVertexSelectionLib RecoToolInterfaces IsolationCorrectionsLib

続いてMyPackageAlg.hでインクルードと宣言を行います。さらに後で使う関数の宣言も行っておきます。

#include "PhotonVertexSelection/IPhotonVertexSelectionTool.h"
#include "PhotonVertexSelection/IPhotonPointingTool.h"
#include "RecoToolInterfaces/ITrackIsolationTool.h"
#include "InDetTrackSelectionTool/IInDetTrackSelectionTool.h"
#include "IsolationCorrections/IIsolationCorrectionTool.h"
...
static bool comparePt(const xAOD::Photon *a, const xAOD::Photon *b);
...
asg::AnaToolHandle<CP::IPhotonVertexSelectionTool> m_vertexTool;
asg::AnaToolHandle<CP::IPhotonPointingTool> m_pointTool;
asg::AnaToolHandle<InDet::IInDetTrackSelectionTool> m_indet_tracksel;
asg::AnaToolHandle<xAOD::ITrackIsolationTool> m_track_iso_tool;
asg::AnaToolHandle<CP::IIsolationCorrectionTool> m_isoCorrTool;
asg::AnaToolHandle<CP::IIsolationSelectionTool> m_photonFCLooseIsoTool;

注意: ミューオンでのisolation選別を行っていない場合は、isolation選別用のライブラリとパッケージのインクルードも行っておきましょう。

MyPackageAlg.cxxでもいくつかインクルードを行います。

#include "xAODEgamma/ElectronxAODHelpers.h"
#include "PhotonVertexSelection/PhotonVertexHelpers.h"

コンストラクタをに以下を追加します。(,を適宜つけるのを忘れないこと)

m_vertexTool("CP::PhotonVertexSelectionTool/vxtTool")
m_pointTool("CP::PhotonPointingTool/pointingTool",this)
m_indet_tracksel("InDet::InDetTrackSelectionTool/trackSel")
m_track_iso_tool( "xAOD::TrackIsolationTool/TrackIsolationTool")
m_isoCorrTool("CP::IsolationCorrectionTool/IsoCorrTool",this)
m_photonFCLooseIsoTool("CP::IsolationSelectionTool/PhotonIsolationSelectionTool", this)

initialize()で初期化を行います。

CHECK( m_vertexTool.retrieve() );
CHECK( m_pointTool.retrieve() );
CHECK( m_indet_tracksel.setProperty("CutLevel" , "Loose") );
CHECK( m_indet_tracksel.setProperty("minPt" , 1000.0 ) );
CHECK( m_indet_tracksel.setProperty("maxZ0SinTheta", 3.0 ) );
CHECK( m_indet_tracksel.retrieve() );
CHECK( m_track_iso_tool.setProperty("TrackSelectionTool", m_indet_tracksel) );
CHECK( m_track_iso_tool.retrieve() );
bool isAtlfast = false;
CHECK( m_isoCorrTool.setProperty( "IsMC", !isData) );
CHECK( m_isoCorrTool.setProperty( "AFII_corr", false) );
CHECK( m_isoCorrTool.retrieve() );
CHECK( m_photonFCLooseIsoTool.setProperty("PhotonWP", "FixedCutLoose") );
CHECK( m_photonFCLooseIsoTool.retrieve() );

さてexecute()でprymary vertex選別を進めていきます。まずはphotonのエネルギー補正や、シャワー補正を済ませておきます。shallow copyに補正済みのphotonを入れておいてください。その上で空のvertexとshallow copyの結果を入れるコンテナを用意します。

const xAOD::Vertex *vertex = nullptr;
xAOD::PhotonContainer signalphotons(photonsSC->begin(), photonsSC->end(), SG::VIEW_ELEMENTS);
signalphotons.sort(comparePt);

ここでphotonsSCはshallow copyのことです。さらにcomparePtは先ほど.hで宣言した関数で、中身は.cxxの一番下に書いておきましょう。

bool MyPackageAlg::comparePt(const xAOD::Photon *a, const xAOD::Photon *b)
{
  return a->pt() > b->pt();
}

コンテナに入っているphotonに対して、preselectionをかけます。どんな選別をかけるかはHGamグループのノートを参照してください。

for (auto photon4vtx = signalphotons.begin(); photon4vtx != signalphotons.end();) {
  if( PASS_PRESELECTION )
    photon4vtx++;
  else
    photon4vtx = signalphotons.erase(photon4vtx);
}

これによってpreselectionをパスしたphotonのみ、コンテナに残ります。
いよいよprimaey vertexを選別します。preselectionをパスしたphotonが2つ未満の場合は、di-photon vertexが探せないので、普通のprimary vertexを選んでおきます。

if(signalphotons.size() < 2) {
  const xAOD::VertexContainer *vertices = nullptr;
  CHECK (evtStore()->retrieve(vertices, "PrimaryVertices"));
  vertex = xAOD::PVHelpers::getHardestVertex(vertices);
}
else {
  signalphotons.resize(2);
  m_vertexTool->getVertex(signalphotons, vertex);
}

選別直前でphotonが2より大きい場合に、強制的に3つ目以降をコンテナから消し去っていることに注意してください。

primary vertexを選び直したので、photonの進行方向 (eta, phi) を計算し直しましょう。これによりdi-photon質量分解能が向上します。

for (auto signalphoton : signalphotons) {
  if (signalphoton->caloCluster() && fabs(signalphoton->caloCluster()->etaBE(1)) < 10.0
      && m_pointTool->correctPrimaryVertex(*signalphoton, vertex->z()).isFailure() ) {
    ATH_MSG_WARNING ("execute(): Pointing tool unable to correct photon for PV");
  }
...

注意: photonの選別において検出器の有効範囲の選別を行いました。しかしその時使用したeta(s2)は単なる検出器上の座標、すなわち検出器中心から見たetaを表しているので、この補正とは無関係です。

さらにこのループ内でisolationの補正も行います。

  std::set<const xAOD::TrackParticle*> tracksToExclude;
  for (unsigned int iv = 0; iv < signalphoton->nVertices(); iv++) {
    const xAOD::Vertex *phvtx = signalphoton->vertex(iv);
    for (unsigned int itk = 0; itk < phvtx->nTrackParticles(); itk++) {
      tracksToExclude.insert(xAOD::EgammaHelpers::getOriginalTrackParticleFromGSF(phvtx->trackParticle(itk)));
    }
  }
  std::vector<xAOD::Iso::IsolationType> m_isoT = {xAOD::Iso::ptcone40, xAOD::Iso::ptcone30, xAOD::Iso::ptcone20};
  xAOD::TrackCorrection m_corrList;
  m_corrList.trackbitset.set(static_cast<unsigned int>(xAOD::Iso::coreTrackPtr));
  m_track_iso_tool->decorateParticle(*signalphoton, m_isoT, m_corrList, vertex, &tracksToExclude);
  if (not m_isoCorrTool->applyCorrection(*signalphoton))
    ATH_MSG_FATAL("execute(): Couldn't correct photon isolation leakage?");

これで補正は終了です。後はイベント選別を行い、di-photonで質量を組めば、Higgsのピークが見つかるはずです。

Overlap Removal

複数のparticleに識別されたobjectのoverlapを解きます。まずCMakeLists.txtを書き換えます。

LINK_LIBRARIES [...] AssociationUtilsLib

続いてMyPackageAlg.hでインクルードと宣言を行います。

#include "AssociationUtils/OverlapRemovalInit.h"
#include "AssociationUtils/ToolBox.h"
...
ORUtils::ToolBox m_orToolbox;

MyPackageAlg.cxxでコンストラクタに以下を追加します。(,を適宜つけること)

m_orToolbox("ORToolbox",this)

initialize()で初期化を行います。今回はb-jetに対して専用のoverlapを使わず、またtauは使用しない設定で行きます。実際は解析に合ったオプションを与える必要があります。

std::string isBTag = ""; // if use Btag OR -> "isBTag" 
ORUtils::ORFlags orFlags("OverlapRemovalTool", "ORInputLabel", "isOR");
orFlags.bJetLabel = isBTag;
orFlags.outputPassValue = false;
orFlags.doElectrons = true;
orFlags.doMuons = true;
orFlags.doJets = true;
orFlags.doTaus = false;
orFlags.doPhotons = true;
orFlags.doFatJets = false;
orFlags.doEleEleOR = false;
orFlags.doMuPFJetOR = false;
CHECK( ORUtils::recommendedTools(orFlags, m_orToolbox));
CHECK( m_orToolbox.muJetORT.setProperty("ApplyRelPt", true) );
CHECK( m_orToolbox.initialize() );

ここで"ORInputLabel"と"isOR"という引数を、orFlagに与えていることを覚えておいてください。
実際に使ってみます。execute()内で各particleのshallow copyを生成し、さらにpreselectionをかけた後で、

CHECK( m_orToolbox.masterTool->removeOverlaps((electrons_shallowCopy.first), (muons_shallowCopy.first), (jets_shallowCopy.first), nullptr, (photons_shallowCopy.first), nullptr) );

のように、ツールにshallow copyを食べさせます。
ここでathenaコマンドを実行してみましょう。残念ながら、エラーを吐いて止まるはずです。

なぜfailしたか、答えは先ほどの"ORInputLabel"にあります。overlap removalは、その解析においてあるparticleらしいと識別されたもののみの間で行われるため、コンテナから取ってきたparticle全ては使用しません。overlap removalにかけるというラベルの役割を果たすのが、"ORInputLabel"です。(別な名前にしたければ、orFlagの引数の文字列を変えればokです。)つまり各particleにラベルを張っていく作業が必要になります。例えばphotonを例にとると:

for (auto photonSC : *photonsSC) {
  ...
  //--- preselection ---
  photonSC->auxdata< char >( "ORInputLabel" ) = pass_preselection;
  ...
}

ここでshallow copyコンテナ内の全てのparticleにラベルの変数を与えること、そしてラベルはcharで与えることに注意してください。
もう一度athenaコマンドを実行してみましょう。今度は成功するはずです。

Overlap removalの結果は、これまたorFlagの引数に入れた"isOR"という名前の変数に格納されます。これはchar型の変数で、removeOverlaps()関数を呼ぶと自動で生成されます。例えば以下のようにすると、結果が確認できます。

for (auto photonSC : *photonsSC) {
  ATH_MSG_INFO("execute(): photon after OR pt/eta/phi/isOR = " << (photonSC->pt() * 0.001) << " GeV/" << photonSC->eta() << "/" << photonSC->phi() << "/" << (int) photonSC->auxdata< char >( "isOR" ));
}

* 注意: muon-jet overlap removalはprimary vertexが組めなかったイベントに対しては正常に動きません。そのようなイベントが気になる場合は、先にprimary vertexがイベントにあることを要求してからoverlap removalを呼びましょう。

MissingET

METを計算します。ここではジェットの pileup rejectionを終えていることを前提として話を進めます。いつもどおり、まずCMakeLists.txtを書き換えます。

LINK_LIBRARIES  [...]  METInterface

続いてMyPackageAlg.hでツールのインクルードと宣言を行います。

#include "METInterface/IMETMaker.h"
...
asg::AnaToolHandle<IMETMaker> m_metMaker;

MyPackageAlg.cxxでもMETコンテナを読めるようにインクルードします。ついでにIParticleHelpersというパッケージもインクルードしておきます。

#include "xAODMissingET/MissingETAuxContainer.h"
#include "xAODMissingET/MissingETAssociationMap.h"
#include "xAODMissingET/MissingETContainer.h"
#include "xAODBase/IParticleHelpers.h"

さらにコンストラクタに以下を追加します。(,を適宜つけること)

m_metMaker("met::METMaker/METMaker")

initialize()でツールの初期化を行います。 ここで"JetRejectionDec"というオプションで指定される変数名は、必ず pileup rejectionの項でfJVTのツールに対して"OutputDec"というオプションで指定した変数名と同一のものにしておきます。

CHECK( m_metMaker.setProperty("ORCaloTaggedMuons", true) );
CHECK( m_metMaker.setProperty("DoSetMuonJetEMScale", true) );
CHECK( m_metMaker.setProperty("DoRemoveMuonJets", true) );
CHECK( m_metMaker.setProperty("JetRejectionDec","passFJvt"));
CHECK( m_metMaker.retrieve() );

execute()ではまずツールを適切に動かすために、setOriginalObjectLink()という関数を呼び出して、shallow copyと元のコンテナの関係を明示しておきます。さらに各オブジェクトのエネルギー等の補正・較正を行った後に、METの計算に使う選別を行います。
選別は解析によって要求を変えますが、一般的にはOverlap Removalの時と同じきつさを要求することが多いです。

bool setLinksElectron = xAOD::setOriginalObjectLink(*electrons,*electronsSC);
if(!setLinksElectron) {
  ATH_MSG_ERROR("execute(): Failed to set original electron links.");
  return StatusCode::FAILURE;
}
for (auto electronSC : *electronsSC) {
  ...
  //--- preselection ---
  electronSC->auxdata< bool >( "pass_baseline" ) = xxx;
}

上の例ではauxdata< bool >( "pass_baseline" )という変数に、選別結果を詰めています。変数名は何でもいいですが、適当な変数にMET用選別をパスしたかどうかの情報を詰めておきましょう。その上でMET計算用のベクターを用意します。

ConstDataVector<xAOD::ElectronContainer> metelectron(SG::VIEW_ELEMENTS);
for (const xAOD::Electron* el : *electronsSC) {
  if( el->auxdata< bool >( "pass_baseline" ) ) metelectron.push_back(el);
}

ミューオン、photon, タウについても同様にしてベクターを用意します。

ジェットについてもやはりsetOriginalObjectLink()を呼び出し、較正等を行っておきます。またJVT, fJVTの計算を済ませておきます。JVTなどの詳しい取り扱いに関してJetETMissグループのrecommendationがあるので、 TWikiを読んでおきましょう。

さてもう少し準備があります。

std::string m_eleTerm = "RefEle";
std::string m_gammaTerm = "RefGamma";
std::string m_tauTerm = "RefTau";
std::string m_jetTerm = "RefJet";
std::string m_muonTerm = "Muons";
std::string m_inputMETMap = "METAssoc_AntiKt4EMTopo";
std::string m_inputMETCore = "MET_Core_AntiKt4EMTopo";
std::string m_outMETTerm = "FinalTrk";
const xAOD::MissingETContainer* metcore(0);
CHECK( evtStore()->retrieve( metcore, m_inputMETCore ) );
const xAOD::MissingETAssociationMap* metMap(0);
CHECK ( evtStore()->retrieve(metMap, m_inputMETMap) );
std::string softTerm = "PVSoftTrk";
metMap->resetObjSelectionFlags();

ここら辺はおまじないだと思ってくれればよいですが、使うジェットの種類やsoft termによっては、書き換えが必要になります。

いよいよMETの計算に取り掛かります。

CHECK( m_metMaker->rebuildMET(m_eleTerm, xAOD::Type::Electron, newMetContainer, metelectron.asDataVector(), metMap) );

以下ミューオン、photon, タウについても同様の手続きを行い、先ほど用意したベクターをMETMakerに食べさせていきます。最後にジェットの計算を行い、METを取り出しましょう。

xAOD::JetContainer metjet(jetsSC->begin(), jetsSC->end(), SG::VIEW_ELEMENTS);
ATH_CHECK( m_metMaker->rebuildJetMET(m_jetTerm, softTerm, newMetContainer, &metjet, metcore, metMap, true) );
m_metMaker->buildMETSum(m_outMETTerm, newMetContainer, (*newMetContainer)[softTerm]->source());

これでnewMetContainerにMETの計算結果が入ります。取り出すときは以下のようにします。

MET = (*newMetContainer)[m_outMETTerm]->met()*0.001;
METx = (*newMetContainer)[m_outMETTerm]->mpx()*0.001;
METy = (*newMetContainer)[m_outMETTerm]->mpy()*0.001;
METphi = (*newMetContainer)[m_outMETTerm]->phi();
METsum = (*newMetContainer)[m_outMETTerm]->sumet()*0.001;

newMetContainerはいらなくなったら消しておきます。

delete newMetContainer;
delete newMetAuxContainer;

Trigger Matching

シミュレーションイベントのevent weighting

MC generator weight

pileup reweight

Run2では、MCサンプルにおけるpileupは digitization (HITS->RDO) の時点でシミュレーションに組み込まれます。しかしMCサンプルを作った時点では、完全にdataのpileup分布を予測できるわけではありません。dataの取得後に予想とのずれを補正するために、イベントのreweightingを行います。ここでは GoodRunsListを使える状態であることを前提としています。

MyPackageAlg.hでインクルードと宣言を行います。

#include "AsgAnalysisInterfaces/IPileupReweightingTool.h"
...
asg::AnaToolHandle<CP::IPileupReweightingTool> m_pileuptool;

MyPackageAlg.cxxのコンストラクタを書き換えます。

MyPackageAlg::MyPackageAlg( const std::string& name, ISvcLocator* pSvcLocator ) :
 AthAnalysisAlgorithm( name, pSvcLocator ),
 ...
 m_pileuptool("CP::PileupReweightingTool/tool", this)

さらにinitilize()で必要なオプションを与えます。

std::vector<std::string> listOfLumicalcFiles = {"GoodRunsLists/data15_13TeV/20170619/PHYS_StandardGRL_All_Good_25ns_276262-284484_OflLumi-13TeV-008.root",
                                                "GoodRunsLists/data16_13TeV/20180129/PHYS_StandardGRL_All_Good_25ns_297730-311481_OflLumi-13TeV-009.root",
                                                "GoodRunsLists/data17_13TeV/20180619/physics_25ns_Triggerno17e33prim.lumicalc.OflLumi-13TeV-010.root",
                                                "GoodRunsLists/data18_13TeV/20190318/ilumicalc_histograms_None_348885-364292_OflLumi-13TeV-010.root"};
std::vector<std::string> configFiles = {"/cvmfs/atlas.cern.ch/repo/sw/database/GroupData/dev/PileupReweighting/mc16_13TeV/pileup_mc16a_dsid343981_FS.root",
                                        "/cvmfs/atlas.cern.ch/repo/sw/database/GroupData/dev/PileupReweighting/mc16_13TeV/pileup_mc16d_dsid343981_FS.root",
                                        "/cvmfs/atlas.cern.ch/repo/sw/database/GroupData/dev/PileupReweighting/mc16_13TeV/pileup_mc16e_dsid343981_FS.root"};
CHECK ( m_pileuptool.setProperty("LumiCalcFiles",listOfLumicalcFiles) );
CHECK ( m_pileuptool.setProperty("ConfigFiles",configFiles) );
CHECK ( m_pileuptool.retrieve() );

ここでconfig fileと呼ばれる、MCのpileupの環境を記録しているファイルは、ggH, H->gamgamサンプルであるmc16_13TeV.343981に対するものを使用しています。他のサンプルを触りたい時は適宜変更してください。

ではexecute()でreweighting factorを計算してみましょう。

float PileupWeight = 1.;
if(isMC) {
  m_pileuptool->apply( *ei );
  PileupWeight = ei->auxdecor<float>("PileupWeight");
  ATH_MSG_INFO("pileup weight=" << PileupWeight );
}

Metadata

サンプルを入手する

rucio

eos

jobをbatchシステムに投げる

login.iceppではbsubを使用できます。使えるキューをbqueuesコマンドで確認してみましょう。

bqueues

実際にキューにjobを投げてみましょう。bsubを付けて実行するだけですが、せっかくなので「テキストファイルにリストされたファイルを一つ一つのjobで回す」スクリプトを書きましょう。例えばこんな内容のスクリプトを、runディレクトリの下に置きます。

#!/bin/bash
BSUB_QUE="4h"
if [ $# != 1 ]; then
  echo "Usage run_bsub.sh <LIST>"
  exit
fi
if [ -f $1 ] ; then
  echo $1
  for FILE in `cat $1`
  do
    TEMP_DIR=`echo "$FILE" | sed -e "s:/:_:g"`
    DIR="Output/$TEMP_DIR"
    mkdir -p $DIR
    cd $DIR
    echo "athena MyPackage/MyPackageAlgJobOptions.py --filesInput=$FILE --evtMax=-1" > tmp_athena_run.sh
    # cat tmp_athena_run.sh 
    bsub -q $BSUB_QUE -o log.out < tmp_athena_run.sh
    cd -
  done
fi

実行してみましょう。

cd $TestArea/../run
ls /gpfs/fs7001/youhei/AnaTutorial/mc16_13TeV.361107.PowhegPythia8EvtGen_AZNLOCTEQ6L1_Zmumu.merge.AOD.e3601_s3126_r9781_r9778/* > list4bsub.txt
./run_bsub.sh list4bsub.txt

走っているjobの様子は

bjobs

あるいはもっと詳しい情報を

bpeek <JOB_NUMBER>

で見ることができます。

jobをGRIDに投げる

以下はGRIDの証明書を持っていることを前提としています。持っていない人は ATLASソフトウェア講習を参照。

GRID上のサンプルに対してjobを走らせたい時は、GRIDにjobを投げるか、GRIDからデータをダウンロードして、ローカルで走らせるかが選択肢として存在します。(例外として、TOKYOサイトに置いたサンプルを、login.iceppのシステムからアクセスする方法もあります。)シグナルサンプルなどはそれほど重くないですし、何度も走らせることが多いので、後者でやることが多いかもしれません。一方実データやBGサンプルは数が多いので、GRIDに投げたくなるでしょう。
GRIDへのjob投入は基本的には簡単です。GRIDに投げるためのpandaと呼ばれるクライアントを立ち上げて、athenaコマンドを少し変えるだけです。

lsetup panda
lsetup rucio
voms-proxy-init -voms atlas
pathena MyPackage/MyPackageAlgJobOptions.py --inDS=INPUT_DATASET --outDS=user.YOURUSERNAME.xxx

(2, 3行目は環境によっては不要です。)INPUT_DATASET, YOURUSERNAMEは適宜書き換えてください。

jobがsubmitされたら、pandaのモニタリングサイトでjobの様子を見てみましょう。
https://bigpanda.cern.ch/tasks/?username=FirstName%20FamilyName
FirstName, FamilyName は忘れずに設定し直してください。
GRID jobは結構ややこしくて、いろいろな理由ですぐ死にます。例えば1 fileにつき1 jobだとjob数が多すぎて怒られることがあります。このような場合は--nFilesPerJobオプションや--nGBPerJobオプションで調整してあげましょう。(例: --nGBPerJob=MAX)

もしかしたらCENTOS7でjobをgridに投げると、いくつかのjobが失敗するかもしれません。そんな時は最初の環境設定で行ったsetupATLASを少し変えて、

setupATLAS -c slc6

として、この環境でjobを投げてみてください。

コマンドライン上でjobを管理することも可能です。失敗したjobをretryしたり、いらないjobをkillしたりできます。 pbookというコマンドを試してみましょう。詳しくは pbookのTWiki参照。

pbook

frameworkをgitで管理

読み物

EventLoop vs Athena

解析手法は歴史的な理由 (後述) により2種類が存在します。

このチュートリアルではAthenaを使用します。最近はほとんど違いが無くなりましたので、あまり気にする必要はありません。具体的にはasetupでの違い、CP toolの導入の違いが主だったものとしてあります。あとはAthenaのほうが使える情報は少し多い、くらいです。
しかしどちらを使っているかがわからないと、TWikiで情報を漁るときに混乱するので注意しましょう。

例えば他のチュートリアルでは、

  • EventLoop: Analysis Software Tutorial
  • Athena: xAOD Analysis in CMakeとかATLASソフトウェア講習
を使用しています。

* 歴史的なお話: 昔xAODがまだAODという名前だったころ、AODはROOTからは読めませんでした。そのためAODを使って解析していた人たちは、Athenaを使用していました。そのころDxAODは存在せず、D3PDというシンプルなnTupleでした。Branchも全てfloatやint、あるいはそのvectorでした。こちらは当然ROOTで読めます。つまり解析の対象のファイルが全く異なるフォーマットだったため、解析frameowrkもかけ離れていたのです。
それがRun2でxAOD/DxAODが導入されて同じような構造になったため、解析frameworkも近しいものになりました。しかし未だ統一はされず、中途半端に分かれています。

ATLASのオフィシャルコードが分からない!

近くのスタッフか先輩に聞くのが有効です。他には、ATLASのTWikiの右上の検索窓から検索をかけて情報を探すのも手です。この時検索先を"Twiki"から"All"にし、さらに"cernscope:Atlas"を消しておくと、過去のCDSのスライドやE-groupsでのやり取りも出てきます。
以下ではそれ以外の方法をいくつか紹介します。

コードを読む

athenaの各パッケージの中身が知りたい場合は、 GitLXRが使えます。

TWikiを読む

CPツールのインストラクションは、各CPグループのTWikiに載っています。各ツールで行っていることや、使い方が説明されています。CPグループのページから辿ってみましょう。

他のframeworkの真似をする

歴史的なお話:昔々Run1の頃、人々は自分でスクラッチから書いたコードを使って、D3PDやAODを解析していました。それでは効率が悪かろうと、いち早く共通のframworkを使い始めたのがSUSYグループです。しかしたびたび起きるAOD/D3PDの仕様変更や、解析ごとに独自のツールを開発していたことなどにより、かなり大変なお仕事でした。その後Run2になると変化が徐々に緩やかになり、解析ごとの手法の違いを減らす努力も働いたため、共通のframeworkを使用するのが一般的になりました。そんなこんなで今日では、多くの解析で似たような解析frameworkを使っています。例えばTopなら AnalysisTop、VH(->bb) なら CxAODというように、それぞれの解析にそれぞれの解析frameworkがあります。困ったことがあれば、それらを参考にするといいでしょう。

メーリスで聞く

専用のE-groupsがいくつか用意されています。例えばATLASのソフトウェアに関することは

https://groups.cern.ch/group/hn-atlas-PATHelp/default.aspx

で聞くと助けてくれます。

ATLASのデータの種類

サンプルはATLAS centralに加工され、解析レベルではxAODまたはDxAODというフォーマットでユーザーに届けられます。

ESDには再構成された全情報が入っています。その下流のxAODではほとんどの解析で使わない情報がそぎ落とされ、軽くなります。この段階からrootで読むことができます。その下流のDxAODではさらにデータは軽くなり、多くの解析ではこのDxAODを使用します。DxAODは解析ごとに用意されており、ユーザーは自分の見たい解析に合ったDxAODを選ぶことが必要になります。またxAOD->DxAODの過程で加えられる情報もあるため、xAODなら全ての解析をカバーできるかというとそうとは限りません。

一般にデータの加工は何度も行われます。検出器の理解が進むことで、より品質の高い、あるいはより実データに近いサンプルを作成するためです。上流ほど生のデータに近いため安定であり、逆に下流に行くほど加工の回数も増えていきます。

実データ

トリガーされたイベントの情報は、RAWデータと呼ばれます。これはそのままでは解析には使えないので、いくつかのステップを経て、解析用のxAOD/DxAODという形式に加工されます。一方RAWデータはテープに記録され、大切に保管されます。

データの加工は繰り返されるので、同じRunのデータがいくつか存在することになります。さらにstreamと呼ばれるトリガーのタイプごとにデータを分けて保存するので、この意味でもデータはいくつか存在することになります。さらにさらに最下流のDxAODは解析の種類によって複数存在するため、結果として1 Runにつき、100を超えるデータが存在します。

データセットの名前の見かた

データセットの名前の例として、チュートリアルで使用するサンプルを挙げます。下はxAODの例。

DxAOD の場合は、オレンジの部分が書き換わり、derivationのtagが末尾に追加されます。

tagの詳細は AMIで調べることができます。

ところでmc16はさらにmc16a, d, e (あとそれ以外) に分かれていて、それぞれ役割が変わります。

campaign 対応するデータ
mc16a 2015 + 2016
mc16d 2017
mc16e 2018
それ以外 Run2のfullの解析では基本的には使わない
これらは tagによって見分けることが可能です。

DxAOD の構造

Monte Carloシミュレーションサンプル

ジェットの種類

実践編のジェットの項目では"AntiKt4EMTopoJets"というジェットを使っています。DxAODには他にも何種類かのジェットが入っています。checkxAOD.pyで確認してみましょう。

checkxAOD.py /gpfs/fs7001/youhei/AnaTutorial/mc16_13TeV.361107.PowhegPythia8EvtGen_AZNLOCTEQ6L1_Zmumu.deriv.DAOD_STDM3.e3601_s3126_r10201_p3526/DAOD_STDM3.13992482._000240.pool.root.1 | grep Jet

track jetやLCTopo jet, PF jetの存在が確認できます。一方xAODにはtrack jetなどは入っていません。これはこれらのジェットがxAODをDxAODに加工する過程で生成されるものだからです。したがってDxAODの種類によって入っているジェットが異なります。

いくつかあるジェットの違いはジェットの再構成アルゴリズムの違いです。例えば"AntiKt4"はΔR 0.4のanti-kTアルゴリズムでジェットが作られたことを意味します。anti-kTアルゴリズムについては各自ググってください。対抗馬にCambridge-Aachenアルゴリズムで作られたCA jetがいます。
それ以外は

  • EMTopoとLCTopoはキャリブレーションの違いです。とりあえずLCTopoには触らないでokです。これらはカロリメータのクラスターによって構成されています。
  • track jetはその名のとおり、trackによって構成されるジェットです。中性粒子の成分は見ていません。
  • PF jetは、particle flowアルゴリズムによって構成されたジェットです。雑にいうと荷電成分はtrackから取ってきて、中性成分をカロリメータから取ってきているハイブリッドなジェットです。
  • truth jetはMC truthの情報から作られたジェットです。検出器が完璧に理想的であれば実現できるジェットとほぼ同義になります。TruthWZと、"WZ"がついているのは、truth jetに対してさらに、通常は落とすWZ由来のミューオン、ニュートリノも使用したtruth jetです。
  • tau jetはタウを捕まえるのに特化したジェットです。詳しくはタウを扱うときに。
実際に解析ではどれを使うのが良いでしょう。
  • 特別なジェットが不要なら、AntiKt4EMTopoJetsでokです。
  • PF jetは今後スタンダードになるかもしれませんが、2019年4月現在まだバグが多いので慌てて使う必要はありません。
  • track jetはジェットの周りに他の粒子がいて分離したい、つまりΔRを下げたい場合に有効です。ただし中性成分を無視するため、energyはunderestimateされます。
  • TopやVectorボソンの崩壊ハドロンを一つのジェットとして捕まえたい(個々に再構成しようにも、中途半端に重なってしまう)ときは、ΔRを大きくして、AntiKt10を使うとよいでしょう。

-- YoheiYamaguchi - 2019-03-26

Topic attachments
I Attachment History Action Size Date Who Comment
PDFpdf 20190412_AnaTutorial_0.pdf r1 manage 681.6 K 2019-04-12 - 04:43 YoheiYamaguchi  
PDFpdf 20190419_AnaTutorial_1.pdf r1 manage 622.9 K 2019-04-19 - 03:11 YoheiYamaguchi  
PDFpdf 20190426_AnaTutorial_2.pdf r2 r1 manage 3766.6 K 2019-06-12 - 04:58 YoheiYamaguchi  
PDFpdf 20190524_AnaTutorial_3.pdf r2 r1 manage 3648.3 K 2019-05-24 - 05:43 YoheiYamaguchi  
PDFpdf 20190531_AnaTutorial_4.pdf r2 r1 manage 975.9 K 2019-05-26 - 15:25 YoheiYamaguchi  
PDFpdf 20190614_AnaTutorial_5.pdf r1 manage 2038.6 K 2019-06-12 - 04:59 YoheiYamaguchi  
PDFpdf 20190621_AnaTutorial_6.pdf r2 r1 manage 1799.9 K 2019-06-21 - 02:17 YoheiYamaguchi  
PDFpdf 20190627_AnaTutorial_7.pdf r1 manage 793.1 K 2019-06-26 - 19:34 YoheiYamaguchi  
PDFpdf 20190719_AnaTutorial_8.pdf r1 manage 1317.8 K 2019-07-17 - 16:52 YoheiYamaguchi  
PDFpdf 20190725_AnaTutorial_9.pdf r1 manage 2842.1 K 2019-07-25 - 05:55 YoheiYamaguchi  
Edit | Attach | Watch | Print version | History: r64 < r63 < r62 < r61 < r60 | Backlinks | Raw View | WYSIWYG | More topic actions
Topic revision: r64 - 2019-07-25 - YoheiYamaguchi
 
    • Cern Search Icon Cern Search
    • TWiki Search Icon TWiki Search
    • Google Search Icon Google Search

    Main All webs login

This site is powered by the TWiki collaboration platform Powered by PerlCopyright & 2008-2019 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding TWiki? Send feedback