ROOT Tutorial per Lab5

Francesco Cirotto - cirotto@naNOSPAMPLEASE.infn.it

Indice

Introduzione

L'idea di questo tutorial è di fornire gli strumenti per eseguire semplici operazioni sugli oggetti fisici di interesse conservati in una nTupla di root. Dovete creare dei tree su cui fare le analisi a partire da file di testo contenenti variabili cinematiche.

In particolare, i file di testo contengono informazioni su un decadimento in due fotoni. Una parte dei decadimenti derivano dal bosone di Higgs, i rimanenti da processi detti di "fondo". Dovrete individuare la corretta distribuzione del fondo ed effettuare un fit alla distribuzione di massa invariante dei vostri eventi.

Settings

Installazione di root

Potete fare riferimento alla guida ufficiale per quanto riguarda l'installazione. Sono previste diverse modalità di installazione a seconda del sistema operativo che si ha a disposizione.

Potrebbe essere utile per il futuro installare conda, date un'occhiata.

Installazione di wsl2 (Opzionale)

Se non lavorate in un ambiente UNIX (Mac o Linux) avrete bisogno di installare un terminale per fare il tutorial. Microsoft mette a disposizione degli utenti di windows 10 l'utility wsl2 ( sottosistema Windows per Linux) che permette di eseguire un ambiente GNU/Linux, inclusi quasi tutti gli strumenti da riga di comando, le utilità e le applicazioni, direttamente in Windows, senza modifiche e senza il sovraccarico di una macchina virtuale tradizionale o di una partizione.

Una guida esaustiva su tutti i passaggi da seguire la trovate qui. In sintesi si deve

• installare wsl2, cioè windows subsystem for linux, che permette di avere un terminale esattamente come in linux
• installare la GUI, cioè l'interfaccia grafica per aprire in maniera grafica i programmi
• installare un X-server, VcXsrv (seguite le indicazioni)

Potete fermarvi qui, senza fare la parte Full Linux Desktop, che permette di installare tutto linux (non solo il terminale) in windows.

Opzionale: Potreste pensare di proseguire con l'installazione di ROOT per lavorare sul vostro pc (o potrebbe tornarvi utile in futuro). Seguite questa guida o questa.

Informazioni preliminari per il tutorial

I file di testo si trovano in allegato alla pagina I file di testo sono composti da 9 colonne, così identificate

mass   Pseudorapidity  Et1     Eta1    Phi1    Et2     Eta2    Phi2    isHiggs

Tranne l'ultima colonna, che è di numeri interi, il resto sono tutti numeri decimali. Nello specifico dunque avrete a disposizione la massa invariante e la pseudorapidità (eta) del sistema dei due fotoni, l'energia trasversa, eta e phi di ciascun fotone della coppia. L'ultima colonna (isHiggs) vi permette di selezionare eventi provenienti dal decadimento dell'Higgs (quando è uguale a 1) da eventi di fondo (quando è 0).

ROOT tips and tricks

• ROOT intro slide
• I tutorial forniti sul sito di ROOT Link ai tutorial
• Breve esempi su come utilizzare la classe degli Istogrammi
• Capitolo sul fit
• Un pò di estetica: legenda, opzioni grafiche, colori, stile
• Classe dei TLorentzVector

Tutorial

Come creare un Tree a partire da un file di testo

Un esempio di come scrivere un file .root contenente un Tree a partire da un file di testo si trova qui. Di seguito una versione molto simile

void writeTree()
{

   //stringhe che rimandano ai nome dei file di output e di input
  TString filenameROOT = "esempio.root";
  TString filenameInput = "ilmioinput.txt"

  TFile* myFile = TFile::Open(filenameROOT,"RECREATE");
  TTree *tree = new TTree("Nome_tree","Qui potete metterci un titolo a scelta");


  tree->ReadFile(filenameInput, "x/D:y/D:z/I");
  tree->Write();
  myFile->Close();

}

Vediamo nel dettaglio alcune righe

TFile* myFile = TFile::Open(filenameROOT,"RECREATE");
  TTree *tree = new TTree("Nome_tree","Qui potete metterci un titolo a scelta");

Queste due righe si occupano di creare un file root, il cui nome è contenuto nella stringa filenameROOT (l'opzione "RECREATE" serve per sovrascrivere ogni volta il file), e di creare un tree con nome e titolo di vostra scelta. La parte di scrittura è gestita da una singola riga

tree->ReadFile(filenameInput, "x/D:y/D:z/I");

dove filenameInput è il file di testo da leggere e il secondo argomento del metodo ReadFile indica il nome (ed il tipo) che attribuite a ciascuna colonna. Nell'esempio dunque, state salvando 3 colonne nel tree con i nomi x,y e z. x e y sono double (D) e z è un int (I). Modificate secondo la formattazione del vostro file questa riga, in base anche alla descrizione fatta precedentemente. Il metodo utilizzerà di default lo spazio come delimitatore tra due colonne. Nota: il metodo ignora nel file di testo tutte le righe che iniziano con il carattere '#'

ATTENZIONE! nel caso in cui avete a che fare con un testo che utilzza un altro delimitatore, ad esempio una virgola, potete specificare questa cosa nel seguente modo

tree->ReadFile(filenameInput, "x/D:y/D:z/I", ',');

Le ultime due righe scrivono l'oggetto tree nel file, e chiudono il file.

Aprire e visualizzare il tree

root -l ilmiofile.root

per vedere cosa c'è all'interno del file root:

.ls

Supponendo che il tree si chiami ilmiotree

per vedere quali variabili sono presenti nel tree:

ilmiotree->Print()

Potete utilizzare il metodo Draw() per visualizzare la distribuzione delle variabili. Il metodo vi permette di applicare anche una selezione. Supponiamo che le variabili salvate abbiano i seguenti nomi mass, Pseudorapidity ,Et1 ,Eta1,Phi1,Et2, Eta2 ,Phi2 ,isHiggs posso graficare la distribuzione della massa facendo

ilmiotree->Draw("mass")

o, con una particolare selezione:

ilmiotree->Draw("mass", "isHiggs==1")
ilmiotree->Draw("mass", "Et1>60 && Et2<45")

Scrivere una macro per creare istogrammi

Per poter riempire degli istogrammi a partire dalle variabili del tree esistono diversi metodi.

Tra i più immediati c'è quello di utilizzare il metodo della classe TTree , Draw salvando il risultato in un istogramma. Vediamo nel dettaglio un esempio (spiegazione riga per riga nei commenti)

void makeHisto(){

  
  TFile *file = new TFile("ilmiofile.root"); //Crea e apre il file root esistente
  TTree *tree = (TTree*)file->Get("ilmiotree"); //Recupera nel file un oggetto di tipo TTree che si chiama ilmiotree

  TH1F *myHisto = new TH1F("h","titolo h", 50, 100,300); //nome e titolo dell'histo, numero di bin, minimo e massimo dell'intervallo
   myHisto->SetLineColor(2); //Colore della linea rosso
   myHisto->SetMinimum(0); //Setta il minimo
   myHisto->GetXaxis()->SetTitle("titolo asse X");
   myHisto->GetYaxis()->SetTitle("titolo asse Y");
 
  tree->Draw("mass>>h","isHiggs==1","goff"); //la distribuzione della variabile mass viene riempita nell'histo con nome h, con la selezione applicata. L'opzione "goff" non visualizza l'istogramma,


  TCanvas *c =  new TCanvas ("c","c", 800,800); //Creo una finestra. Se non lo faccio quando utilizzo il metodo Draw dell'histo viene comunque creata di default
  myHisto->Draw();

  c->Print("ilmioplot.pdf"); //salvataggio in pdf dell'istogramma
  myFile->Close();

}

ATTENZIONE! questo metodo funziona molto bene, ed è veloce nel caso in cui le selezioni da applicare sono semplici, e riconducibili alle variabili salvate nel tree.

Un secondo metodo per analizzare i dati è quello di accedere alle variabili evento per evento e scrivere la propria selezione.

void tree1r(){
   // Supponiamo di star leggendo un file col nome "tree1.root", contenente un tree col nome "t1"

   TFile *f = new TFile("tree1.root");
   TTree *t1 = (TTree*)f->Get("t1");

   Float_t px, py, pz; //Devo creare le variabili che punteranno a quelle salvate nel tree

   //Punto le variabili create precedentemente ai branch del mio tree
   t1->SetBranchAddress("px",&px); //quindi la variabile px punta al branch con nome "px"
   t1->SetBranchAddress("py",&py);
   t1->SetBranchAddress("pz",&pz);

   // creo un histo
   TH1F *hpx   = new TH1F("hpx","px distribution",100,-3,3);

   //Accedo evento per evento alle variabili ed eseguo la mia analisi

   Int_t nentries = (Int_t)t1->GetEntries(); //Questo metodo mi restistuisce il numero di eventi salvati nel tree

  //Eseguo un loop
   for (Int_t i=0; i<nentries; i++) {
      t1->GetEntry(i); //questo metodo "carica" l'evento i-esimo con tutti valori dei branch

     //posso ora scrivere la mia selezione
      if(px>0 && py<1){
      hpx->Fill(px); //Utilizzo il metodo Fill per riempire l'istogramma 
      }
   } //Fine loop

   hpx->Draw();  
}

Un'alternativa all'esempio di sopra è quello di utilizzare il metodo MakeClass. Questo genera una class personalizzata che permette di accedere ad un loop pre-scritto da ROOT, dove implementare l'analisi. Può essere utile se si ha a che fare con tree con molte variabili, dove il SetBranchAddress va fatto molte volte. Per maggiori informazioni potete guardare nella guida di ROOT .

Eseguire un fit

Vogliamo ora eseguire un fit ad un istogramma. Anche qui si può proseguire in modi diversi a seconda dei casi. Supponiamo che myHisto sia l'istogramma sul quale vogliamo eseguire un fit.

Vediamo alcuni esempi:

TF1 *gaussian = new TF1("gaussian", "gaus", 80, 100);  //Definizione di una funzione gaussiana nell'intervallo 80,100. 
myHisto->Fit("gaussian","R");

TF1 *pol2 = new TF1("parabola", "pol2", 80, 100);  //Definizione di una parabola nell'intervallo 80,100. 
myHisto->Fit("parabola","R");

TF1 *pol2_V2 = new TF1("parabolaV2", "[0] + [1]*x + [2]*x*x", 80, 100);  //Definizione di una parabola nell'intervallo 80,100. 
myHisto->Fit("parabolaV2","R");

ROOT fornisce, per le funzioni più comuni, delle funzioni predefinite. Negli esempi precedenti ; gaus e po2 (in generale polN, dove N va specificato). Si veda qui per una lista completa. Nel terzo esempio, la funzione è stata scritta utilizzando una formula, con i parametri [0],[1],[2] da stimare nel fit (si noti che è sempre una parabola). L'intervallo di definizione nell'esempio è 80,100. L'opzione "R" nel metodo "Fit" sta ad indicare di eseguire il fit nel range di definizione della funzione. Altre opzioni per il fit qui.

Può essere utile "aiutare" il fit settando opportunamente i parametri. Con riferimento all'esempio della gaussiana, abbiamo 3 parametri da stimare, nell'ordine costante di normalizzazione, media e sigma. Possiamo pensare di settare i parametri nel seguente modo

gaussian->SetParameters(100,5,1); //Stiamo settando nell'ordine costante, media e sigma

oppure pensare di agire sul singolo parametro:

gaussian->SetParameter(1, 125); //Setta il parametro 1 della gaussiana (la media) a 125
gaussian->SetParLimits(2, 2,10); // Limita il parametro 2 nell'intervallo 2,10

Possiamo pensare di fissare un parametro del fit, e quindi eliminarlo dalla stima

gaussian->FixParameter(1, 125); //Fissa a 125 il parametro 1

Si può accedere ad un parametro nel seguente modo

Double_t media = gaussian->GetParameter(1); 
Double_t errore_media = gaussian->GetParError(1); 

Double_t integral = gaussian->Integral(20,40)/myHisto->GetBinWidth(1); //Integrale nell'intervallo 20,40. Per renderlo indipendente dall'istogramma va diviso per la larghezza del bin

Possiamo pensare di sommare due funzioni nel seguente modo

TF1 *f_somma= new TF1("f_somma", "[0] + [1]*x + [2]*x*x + gaus(3)", 80, 100);  //Notare che la gaussiana parte dal parametro 3, quindi avrà i parametri 3,4,5
TF1 *f_sommaV2= new TF1("f_sommaV2", " gaus(0) + pol1(3) + gaus(5)", 80, 100);  //Tra parentesi è indicata la prima posizione del parametro della funzione

Esercizi

• Scrivere un tree a partire dal file di testo
• Contare il numero di eventi di segnale, e di fondo dal tree

• Riprodurre la distribuzione della massa invariante per il segnale, per il fondo e per il totale dei due

• Effettuare un fit alla distribuzione della massa invariante per il fondo, segnale e totale. Che funzione avete utilizzato?. Stimare il numero di eventi

//Sia myHisto l'histo totale
//supponiamo che la funzione totale sia la somma di una gaussiana e di un polinomio di primo grado

TF1 * f_totale = new TF1("f_totale","gaus(0) + pol1(3)", 50,200);
//Setting di parametri

myHisto->Fit("f_totale","R");

//Recupero i valori dei parametri per la pol1
Double_t p0 = f_totale->GetParameter(3);
Double_t p1 = f_totale->GetParameter(4);

//Creo una nuova funzione pol1 con i pametri fittati, cambio il colore e la disegno nella stessa canvas
TF1 *fondo = new TF1("fondo", "pol1",50,200);
fondo->SetParameters(p0, p1); //Settati i prametri
fondo->SetLineColor(4); //cambio colore

fondo->Draw("same"); //l'opzione same disegna nella stessa canvas

Alcune info sulla classe vector di C++

#include <vector>
int main()
{
 // Un vettore dinamico
 std::vector <double> v_dyna;

 // Un vettore con 10 elementi (inizializzati tutti a 0)
 std::vector <int> vint_10_elements(10);

 // Un vettore con 10 elementi, inizializzati a 90
 std::vector <int> vint_10(10, 90);

 // Un vettore inizializzato con il contenuto di un altro
 std::vector<int> vcopy(vint_10_elements);

 // Un vettore con 5 valori presi da un altro
 std::vector<int> vcopy_5(vint_10.begin(),vint_10.begin()+5);

}

Se voglio inizializzarlo con dei valori iniziali, posso scrivere così:

std::vector<int> v = {1, 2, 3, 4, 5};
//(attenzione: questo tipo di inizializzazione è disponibile solo nel nuovo standard C++11).

Per accedere all’i-esimo elemento, si può fare in diversi modi:

std::cout << v[2] << std::endl; //come negli array
std::cout << v.at(2) << std::endl;

Anche con il vettore, se cercate di accedere un elemento al di fuori della dimensione dell’array, otterrete un errore.

Per aggiungere un elemento ad un vettore basta fare

v.push_back(10);

// per avere la dimensione
std::cout << v.size() << std::endl;

Un loop con i vettori si può fare in diversi modi:

std::vector<int> v = {0, 1, 2, 3, 4, 5};

// ciclo for con iteratori
for (std::vector<int>::iterator it = v.begin(); it != v.end(); ++it)
{
std::cout << *it << "\n"; //notate che il contenuto del vettore è *it
}


 // iterazione basata su indici
for (int i = 0; i<v.size();i++)
{
std::cout << v.at(i) << "\n"; //si accede all'i-esimo evento
}

Altre info sull'utilizzo dei vector in c++ qui