Counting mass of TTbar from Monte Carlo generator

Data preparing

First of all we need is make data wich we could operate with. We will use pythia.

Generating pp simulation for TTbar

Athena jobOption setting.

csc_evgen_trf.py [options] <runnumber> <firstevent> <maxevents> <randomseed> <jobconfig> <outputevgenfile> [histogramfile] [ntuplefile] [inputgeneratorfile]

Arguments:
    1 runNumber (int) # each run number corresponds to one physics process
    2 firstEvent (int) # number of the first event in the output data file
    3 maxEvents (int) # Maximum number of events to process
    4 randomSeed (int) # random seed for physics generators
    5 jobConfig (list) # jobOptions fragment containing the physics and the configuration settings
    6 outputEvgenFile (str) # Output file that contains generated events
  [ 7 histogramFile] (str) default='NONE' # Output file that contains histograms.
  [ 8 ntupleFile] (str) default='NONE' # Output file that contains ntuples.
  [ 9 inputGeneratorFile] (str) default='NONE' # Input file used by the particle generator to generate events

$ source setup.sh -tag=14.1.0.3,32,setup,AtlasProduction,releases
$ csc_evgen_trf.py 1 0 10000 1 

TODO

Making .aan.root file from .pool.root

athena -trf

TODO

Data procesing

Making script in root

Now we could make some calculations with data on generator level. For example we could count mass of W from W->lnu decay. We take file from previous chapter .aan.root and open it in root. Go to directory where you have file "TTbar100k.aan.root" and run root. In root open root browser. For example:

root [0] TBrowser b

Warning, important Note that you need to have configured athena (lxplus) or installed root (local).

Now in folder tree click on folder where you have the file and double-click on that file. Now in folder tree click on the folder ROOT Files. Here you see all root files wich has been loaded. Choose your .aan.root file, and CollectionTree folder. Now you see branches of event tree. In each branch is one specific information about all particles in all events. For example "PartPx" is branch wich contain x component of momentum vector of all particles. So now we make script for work with these data.

Leave root with command .q and make new .C file (e.g. FindWMass.C ), and write this into it:

{
TFile* file = new TFile("TTbar100k.aan.root");            //open file containing root trees
TTree* tree = (TTree*)file->Get("CollectionTree");        //get tree named "CollectionTree" from file

int nEvts = tree->GetEntries();
cout << "nEvts getted" <<endl;                            //has to be here, otherwise root dies in horror, no idea why
vector<double> *px=NULL,*py=NULL,*pz=NULL,*E=NULL;        //vectors for momenta, etc.
vector<int> *PDG=NULL;
TBranch *bpx=0,*bpy=0,*bpz=0,*bE=0,*bPDG=0;

tree->SetBranchAddress("PartPx",&px,&bpx);                //select branches
tree->SetBranchAddress("PartPy",&py,&bpy);
tree->SetBranchAddress("PartPz",&pz,&bpz);
tree->SetBranchAddress("PartE",&E,&bE);
tree->SetBranchAddress("PartType",&PDG,&bPDG);

TLorentzVector pVector;                                  //Lorentz vector of particle

TCanvas *c1 = new TCanvas("Histogram","Histogram");      //canvas for histogram
TH2F *partEtaPhi = new TH2F("partEtaPhiHist","Histogram of particles", 50, -10, 10, 50, -4, 4); //2D eta,phi histogram of particles
partEtaPhi->SetXTitle("#eta");
partEtaPhi->SetYTitle("#phi"); 




for (int i=0;i<10;i++) {                                                                 //loop over first 10 events
  tree->GetEntry(i);                                                                     //select i event
  
  partEtaPhi->Reset();                                                                        //reset histogram

  cout << "=========================== " << i << " ==================================="<<endl;

  for (int j = 0; j < px->size(); ++j) {                                                   //loop over all particles in event
    pVector.SetPxPyPzE(px->at(j),py->at(j),pz->at(j),E->at(j));
    partEtaPhi->Fill(pVector.Eta(),pVector.Phi(),pVector.Perp()); 
  }
  partEtaPhi->Draw("LEGO1");
  c1->Update();
  getchar();
}

}

Run root with command

root -x FindWMass.C

Root make new window with 2D histogram of eta, phi with pt weight from all particles in one event. The highest bricks have highest transverse momentum. For more information about eta, phi and pt visit site. Go back to root console and press enter to move to the next event. You could see that historgram is changing. After 10th event procedure ends. If this works, you could now quit root and make some counting. Change your file FindWMass.C to this

{
gROOT->Reset();
#include "TLorentzVector.h"
#include <vector>

//===================================================================================
// OPENING FILE AND TREE
//===================================================================================
TFile* file = new TFile("TTbar100k.aan.root");            //open file containint root trees 
TTree* tree = (TTree*)file->Get("CollectionTree");        //get tree named "CollectionTree" from file 

int nEvts = tree->GetEntries();
cout << "nEvts getted" <<endl;                            //has to be here, otherwise root dies in horror
vector<double> *px=NULL,*py=NULL,*pz=NULL,*E=NULL;        //vectors for momenta, etc.
vector<int> *PDG=NULL;                                    //vectors for PDG
TBranch *bpx=0,*bpy=0,*bpz=0,*bE=0,*bPDG=0;               //branch decalrations

tree->SetBranchAddress("PartPx",&px,&bpx);                //select branches
tree->SetBranchAddress("PartPy",&py,&bpy);
tree->SetBranchAddress("PartPz",&pz,&bpz); 
tree->SetBranchAddress("PartE",&E,&bE);
tree->SetBranchAddress("PartType",&PDG,&bPDG);




//===================================================================================
// VARIABLES AND HISTOGRAM
//===================================================================================
TLorentzVector pVector, v1,v2;               // Lorentz vectors of particle (pVector), leptons (lVec), neutrinos (nuVec) and help vectors (v1,v2)
vector<TLorentzVector> lVec;    
vector<TLorentzVector> nuVec;

for(int k=0; k<4; k++){lVec.push_back(TLorentzVector(0.,0.,0.,0.)); nuVec.push_back(TLorentzVector(0.,0.,0.,0.));} //set four zero Lorentz vectors for leptons and neutrinos 

TCanvas *c1 = new TCanvas("Histogram","Histogram");
TH1F *WMassHist = new TH1F("WMassHist","Histogram of W mass", 200, 35e3, 145e3); 
WMassHist->SetXTitle("mass [MeV]");
WMassHist->SetYTitle("count");




//=================================================================================== 
// LOOPING EVENTS AND PARTICLES
//===================================================================================
for (int i=0;i<10000;i++) {                                                                 //loop over all events
  tree->GetEntry(i);                                                                     //select i event 

  cout << "=========================== " << i << " ==================================="<<endl;
  pVector.SetPxPyPzE(0.,0.,0.,0.);
  for(int k=0; k<4; k++){lVec[k].SetPxPyPzE(0.,0.,0.,0.); nuVec[k].SetPxPyPzE(0.,0.,0.,0.); /*cout<<lVec[k].Mag()<<nuVec[k].Mag()<<endl;*/} 

  for (int j = 0; j < px->size(); ++j) {                                                   //loop over all particles in event
    pVector.SetPxPyPzE(px->at(j),py->at(j),pz->at(j),E->at(j));

    //filter
    if ( (pVector.Perp()  > 150. ) &&                                    //pt cut 
        ((TMath::Abs(PDG->at(j))==13)||                                  //muon/antimuon filter
    (TMath::Abs(PDG->at(j))==14)||                                  //muon neutrino/antineutrino filter
    (TMath::Abs(PDG->at(j))==11)||                                  //electron/positron filter 
    (TMath::Abs(PDG->at(j))==12) )  ){                              //electron neutrino/antineutrino filter

     // lepton-neutrino pair detection
     switch(PDG->at(j)){                                             
       case  13: if (lVec[0].Mag()==0)  { lVec[0] = pVector; break;} //muon
       case -14: if (nuVec[0].Mag()==0) { nuVec[0]= pVector; break;} //muon antineutrino

       case -13: if (lVec[1].Mag()==0)  {lVec[1] = pVector; break;} //antimuon
       case  14: if (nuVec[1].Mag()==0) {nuVec[1]= pVector; break;} //muon neutrino 

       case  11: if (lVec[2].Mag()==0)  {lVec[2] = pVector; break;} //electron
       case -12: if (nuVec[2].Mag()==0) {nuVec[2]= pVector; break;} //positron neutrino

       case -11: if (lVec[3].Mag()==0)  {lVec[3] = pVector; break;} //positron
       case  12: if (nuVec[3].Mag()==0) {nuVec[3]= pVector; break;} //electron neutrino 
       default: break;
     }
   } 

  } 



//===================================================================================
// W MASS COUNTING
//===================================================================================
  for (int k=0; k<4; k++){                                                      
    v1= lVec[k]; v2 = nuVec[k];                               //save pair to help Lorentz vectors for easy to counting
    if ( v1.Mag() * v2.Mag() != 0.) {                         //is needed to make counting?
      cout << ( v1 + v2 ).Mag() << endl;                      //write W mas
      WMassHist->Fill((v1+v2).Mag());                         //add mass to histogram
    }
  }
}



//===================================================================================
// PLOT THE HISTOGRAM
//===================================================================================
WMassHist->Draw();
c1->Update();
}

and run root with same command as before. On the screen you will see mass of W in keV every event every event separately. For more information about used function, visit sites about TLorentzVector, TTree, vector.

Makefile and compiling

Next calculation will need more complex solution. First of all we need to install SpartyJet (manual) tool for root. Then copy all header files (".h") into your include folder (e.g. /usr/include/). Make new ".C" file (e.g. FindTopMass) in directory with ".aan.root" fiel and write write this code:

findTopMass.C

We need make makefile.

Warning, important Warning! Some editors have diferent ! interpretation. Edit this file only in special text procesors (e.g. mcedit).

TODO

Reconstruction and ATLAS simulation

TODO


-- JakubCuth - 28 Sep 2009

Topic attachments
I Attachment History Action Size Date Who Comment
Unknown file formatgz spartyjet_v2.20.0.tar.gz r1 manage 4253.2 K 2009-08-03 - 21:52 JakubCuth This is version wich i have used. Search the web for new version.
Edit | Attach | Watch | Print version | History: r7 < r6 < r5 < r4 < r3 | Backlinks | Raw View | WYSIWYG | More topic actions
Topic revision: r7 - 2009-09-28 - JakubCuth
 
    • Cern Search Icon Cern Search
    • TWiki Search Icon TWiki Search
    • Google Search Icon Google Search

    Sandbox All webs login

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