基本的にはROOTができれば大丈夫。

ATLASのデータセット(dataとmc)はTTreeでrootfileにつまっている、そのTTreeをC++等のコードで読んでヒストグラムを作って、絵を作る。

最近の解析はD3PDというフォーマット(中身がTTreeのrootfile)で行うことを前提にいろいろなツールが作られているので、TTreeを使えるようになってください。

このチュートリアルはlogin.iceppのSL6のマシンで行うことを前提としてます。

適当なディレクトリを作ってそのの下で作業をしましょう。

Higgs->WW->lnuqq

ここではHiggsがWWに崩壊して、一つのWがレプトニックに崩壊するモードのサンプルを使ってお勉強をしましょう。 サンプルは

/gpfs/fs3001/atljphys/datafiles/NTUP_COMMON/mc12_8TeV/mc12_8TeV.167930.PowhegPythia8_AU2CT10_ggH400CPS_WWlepnuqq.merge.NTUP_COMMON.e1617_s1581_s1586_r3658_r3549_p1575/

にあります。この中の適当なファイルを一つROOTで読んで中のTTreeを読むクラスを作りましょう。

 $ root -l /gpfs/fs3001/atljphys/datafiles/NTUP_COMMON/mc12_8TeV/mc12_8TeV.167930.PowhegPythia8_AU2CT10_ggH400CPS_WWlepnuqq.merge.NTUP_COMMON.e1617_s1581_s1586_r3658_r3549_p1575/NTUP_COMMON.01343854._000028.root.1

Applying ATLAS style settings...

root [0]
Attaching file NTUP_COMMON.01343854._000028.root.1 as _file0...
Warning in <TClass::TClass>: no dictionary for class AttributeListLayout is available
Warning in <TClass::TClass>: no dictionary for class pair<string,string> is available
root [1] .ls
TFile**         NTUP_COMMON.01343854._000028.root.1     StreamNTUP_COMMON
 TFile*         NTUP_COMMON.01343854._000028.root.1     StreamNTUP_COMMON
  KEY: AttributeListLayout      Schema;1
  KEY: TDirectoryFile   physicsMeta;1   physicsMeta
  KEY: TDirectoryFile   Lumi;1  Lumi
  KEY: TTree    CollectionTree;1        CollectionTree
  KEY: TTree    physics;1       physics
  KEY: TTree    physicsTrigDec;1        physicsTrigDec
root [2] physics->MakeClass();
Info in <TTreePlayer::MakeClass>: Files: physics.h and physics.C generated from TTree: physics
root [3] .q

physics.hとphysics.Cというファイルができているはずです、中身を見てみてください。

たくさん、ブランチが定義されています、ブランチの名前を見るといくつかの種類に分けられます、

まあ、実際の解析で普通に使うのはこれくらいです。

このphysicsというクラスを使ってイベント毎に解析をします。

MyAna

Step1

ここからはC++のコードを書いてきます。

physicsを継承させたMyAnaというクラスを作ります。

ここに書くのは面倒なので/home/morinaga/Tutorial2014/step1においたのでコピーしてください。

Makefile見ればわかりますが、physicsクラスは巨大なのでコンパイルは一度だけ初めにやって、共有ライブラリ作ります。

使い方は、

 $ make physics
 Compiling physics/physics.o 
Compiling Shared physics 
 $ make clean;make
Now Cleaning
 Compiling main/main.o 
 Compiling MyAna/MyAna.o 
 Compiling all Succesfully Done!!!

とすると、anaという人が出来上がるので、それを

./ana -i input.txt

とするとたくさん数字が出てきます。MyAnaの中身をみればわかりますが（わからなかったらやばいよやばいよ）、electronの数を表示しているだけです。 （コンパイラはclang++にしていますが適宜好きなものにしてください、好みです。Makefile変えてね。） ここまでがStep1です。お疲れ様です。

Step2

Electronの数を表示するだけではだめなので、いろいろやりましょう。

ATLASでは、各物理オブジェクトを再構成してます（ジェットとかミューオンとか）、これには偽物もかなり入っているので（step1は信号なのに電子がたくさんいましたね）こいつらを落としていきます。

Electronの選択

電子は偽物多いです、Fake electronとか呼びますがそのほとんどがジェット（パイオン）由来です。これらを落とすためにATLASではいろいろやっています。

step2といいうディレクトリを見てください、さっきのMyAnaのExecute()にElectronのスタンダードなカットを書いておきました。

  // Electron Selection
  std::vector<TLorentzVector> GoodElectrons;GoodElectrons.clear();
  for (int i = 0; i < el_n; ++i) {
    TLorentzVector this_el = TLorentzVector();
    this_el.TLorentzVector::SetPxPyPzE(el_px->at(i)/1000.,el_py->at(i)/1000.,el_pz->at(i)/1000.,el_E->at(i)/1000.);
    if ((el_author->at(i) == 1) || (el_author->at(i) == 3)) {
      if (fabs(el_charge->at(i)) == 1) {
        if ((el_OQ->at(i)&1446) == 0) {
          if(fabs(this_el.Eta()) < 1.37 || (fabs(this_el.Eta()) > 1.52 && fabs(this_el.Eta()) < 2.47)){
            if(fabs(el_trackz0pvunbiased->at(i)*sin(el_tracktheta->at(i))) < 0.5){
              if(fabs(el_trackd0pvunbiased->at(i)/el_tracksigd0pvunbiased->at(i)) < 6.0){
                if(this_el.Pt() > 25.){
                  if(el_tightPP->at(i) == 1){
                    if( el_ptcone20->at(i)/1000./this_el.Et() < 0.14){
                      if( el_Etcone20->at(i)/1000./this_el.Et() < 0.15){
                        GoodElectrons.emplace_back(this_el);
                      }
                    }
                  }
                }
              }
            }
          }
        }
      }
    }
  }
  std::cout << "#of Electron Raw/Select is " << el_n << "/" << GoodElectrons.size() << std::endl;

上から説明してくと、

まず、セレクションを通ったelectronの箱をvectorで作っておきます。ここで、TLorentzVectorというクラスを使っていますが、まあ知ってますよね、しらないなら調べてね。

再構成された電子のすべてをループでみます。まずTLorentzVectorに値をつめます。（ATLASのエネルギー等の標準の単位はMeVなので1/1000してGeVにします。解析に慣れたいまでもたまに忘れるので覚えておいてください）

ここから合計10つのカットをapplyします、

すべてのセレクションを通ったものをさっきの箱に詰めていきます。（emplace_backが使えない人はpush_backを使いましょう。）

最後にいくつの電子がセレクション通ったかを表示します。結構、落ちていることがわかります。

Muonの選択

続いてMuonのセレクションです。

  // Muon Selection
  std::vector<TLorentzVector> GoodMuons;GoodMuons.clear();
  for ( int i = 0; i < mu_staco_n; ++i) {
    TLorentzVector this_mu = TLorentzVector();
    this_mu.TLorentzVector::SetPxPyPzE(mu_staco_px->at(i)/1000.,mu_staco_py->at(i)/1000.,mu_staco_pz->at(i)/1000.,mu_staco_E->at(i)/1000.);
    if(mu_staco_author->at(i) == 6){
      if(fabs(mu_staco_eta->at(i)) < 2.5){
        if(mu_staco_isCombinedMuon->at(i) == 1) {
          if (mu_staco_expectBLayerHit->at(i) == 0.0 || mu_staco_nBLHits->at(i) > 0) {
            if (mu_staco_nPixHits->at(i) + mu_staco_nPixelDeadSensors->at(i) > 0) {
              if (mu_staco_nSCTHits->at(i) + mu_staco_nSCTDeadSensors->at(i) > 4) {
                if (mu_staco_nPixHoles->at(i) + mu_staco_nSCTHoles->at(i) < 3) {
                  int MuHitN = mu_staco_nTRTHits->at(i) + mu_staco_nTRTOutliers->at(i);
                  bool isOKTRT = false;
                  if ((fabs(mu_staco_eta->at(i)) > 0.1) && (fabs(mu_staco_eta->at(i)) < 1.9)){
                    if((MuHitN > 5) && (mu_staco_nTRTOutliers->at(i) < 0.9*MuHitN)){ isOKTRT = true; }
                  }
                  if((fabs(mu_staco_eta->at(i)) <= 0.1) || (fabs(mu_staco_eta->at(i)) >= 1.9)){
                    if(MuHitN > 5){
                      if(mu_staco_nTRTOutliers->at(i) < 0.9*MuHitN){ isOKTRT = true; }
                    }else{
                      if(MuHitN > 5){
                        if(mu_staco_nTRTOutliers->at(i) < 0.9*MuHitN){isOKTRT = true; }
                      }else{ isOKTRT = true; }
                    }
                  }
                  if(isOKTRT){
                    if(this_mu.Pt() > 25.){
                      if (fabs(mu_staco_trackz0pvunbiased->at(i)*sin(mu_staco_tracktheta->at(i))) < 0.5){
                        if(fabs(mu_staco_trackd0pvunbiased->at(i)/mu_staco_tracksigd0pvunbiased->at(i)) < 3.5){
                          if(fabs((mu_staco_id_qoverp->at(i)-mu_staco_me_qoverp->at(i))/(sqrt(mu_staco_id_cov_qoverp_exPV->at(i)+mu_staco_me_cov_qoverp_exPV->at(i)))) < 5.0){
                            if( mu_staco_ptcone20->at(i)/1000./this_mu.Pt() < 0.14){
                              if( mu_staco_etcone20->at(i)/1000./this_mu.Pt() < 0.15){
                                GoodMuons.emplace_back(this_mu);
                              }
                            }
                          }
                        }
                      }
                    }
                  }
                }
              }
            }
          }
        }
      }
    }
  }
  std::cout << "#of Muon Raw/Select is " << mu_staco_n << "/" << GoodMuons.size() << std::endl;

最後に通ったものを詰めて数を表示します。

ジェットの選択

ジェットです、ATLASでは標準的にAntiKtというアルゴリズムでジェットを再構成しています。4という数字はdR < 0.4の大きさを表します。LCTopoというのはざっくりいうとエネルギーの測り方です。気になれば調べてね。

  //Jet Selection
  std::vector<TLorentzVector> GoodJets;GoodJets.clear();
  for (int i = 0; i< jet_AntiKt4LCTopo_n; i++) {
    TLorentzVector this_jet = TLorentzVector();
    this_jet.TLorentzVector::SetPtEtaPhiE(jet_AntiKt4LCTopo_pt->at(i)/1000.,jet_AntiKt4LCTopo_eta->at(i),jet_AntiKt4LCTopo_phi->at(i),jet_AntiKt4LCTopo_E->at(i)/1000.);
    if(this_jet.Pt() > 30.){
      float const_eta = jet_AntiKt4LCTopo_constscale_eta->at(i);
      if (fabs(const_eta) < 4.5) {
        if (this_jet.Pt() < 50. && fabs(const_eta) < 2.4 && fabs(jet_AntiKt4LCTopo_jvtxf->at(i)) < 0.5) {
          if ((jet_AntiKt4LCTopo_isUgly->at(i) == 0) && (jet_AntiKt4LCTopo_isBadLooseMinus->at(i) == 0) ){
            GoodJets.emplace_back(this_jet);
          }
        }
      }
    } 
  }  
  std::cout << "#of Jet Raw/Select is " << jet_AntiKt4LCTopo_n << "/" << GoodJets.size() << std::endl; 

MissingETの再構成(METだよ)

METの再構成は、概念的には簡単です。すべての物理オブジェクトのベクトル和の逆を取ればよいのです。

言葉では簡単でも、実際に行うのはとっっっっても難しいです。METグループがある程度計算してくれますが、解析によって少し修正しなければいけません。

それは難しいので、ここでは示しませんが解析やってればいずれはやります。

  // MET Reconstruction
  TVector2 MET;
  MET.Set(MET_RefFinal_STVF_etx,MET_RefFinal_STVF_ety);

最後にコンパイルして、走らせてどれくらい数が変わるか見てみましょう。

step2はおしまいです。

Step3

ここではstep2でつくったElectron/Muon/Jet/METを使って実際に分布を作ってみましょう。

lnuqqなので、見えるものとしてはレプトン(el or mu)が1つ、ジェットが２本とMETです。

step2ではのべませんでしたが、これらを選ぶときに重要なカットがまだあります、Overlap removalと僕達は呼んでいますが、それぞれの物理オブジェクトが重なっていないということをチェックします。

step3のディレクトリをコピーとかしましょう。コードが新しくなっています。

Overlap Removal

物理オブジェクトのオーバーラップを解くことは大事です、実は電子なものがジェットとして再構成されていた、というのはよくあります。

基本的にはひとつひとつ他のオブジェクトとのdR が近くないかを調べます。

ATLASの解析でよくあるのは、Muon --> Electron --> Jet の順番に、dR < 0.2になっていないかをチェックします。これを簡単に実装しておいたのでMyAna見てください。 （C++11が使えない人は自分で書きなおしてね）

これをおこなったあとに分布をチェックしていきます。

(ヒストグラムを)つくってあそぼ

ヒストグラムくらいは作ったことあると思うので詳しくは述べませんが、いろんな分布を作ってください。例えば、レプトンのPt/Eta/Phiとか、ジェットのPt/Eta/Phiとか、METとか。。。

WW->lnuqqで重要な変数を幾つか挙げておくので、自分で考えて絵を作って見てください。

・W- > lep + nu の 横方向質量(transverse mass)とPt

・W -> qq(jet + jet) の 質量とPt

・WWシステムの質量とかPtとか。（W massを仮定して2次方程式を解きましょう）

この辺の絵を作って誰かに見せましょう。大丈夫そうだったら前のページのDataAndMc見てください。