Taller 3

Convenciones utilizadas en esta guía

• Comandos escritos en la Linea de Comando de Linux (CLI) van en rojo y en marcados

• Texto en azul es para ser reemplazado por algún un valor indicado
• Argumentos opcionales o pasos opcionales van en rosado.

• Lineas de codigo de algun archivo mas en este fondo beige

• Usaré como editor gedit. Si ustedes son familiares con algun otro, pueden usarlo.

Material

El material para este taller, se puede descargar del siguiente link:

• ejercicios-Parte3.tgz: Ejercicios Parte TC3

tar xzvf ejercicios-Parte3.tgz
cd ejercicios-Parte3

En esta carpeta, se encuentra el código necesario para la hacer los ejercicios. También existe una carpeta solución que contiene el código completo.

Ejercicio 1 - Creacion de una clase para hacer el analisis de datos

En este ejercicio daremos un paso atras: hemos visto que los datos de la simulación en Pythia se guardan en un archivo de salida de extensión .root. Estos datos se encuentran en formato especial, un arbol o TTree de ROOT. El análisis se realizo mediante una clase, la cual es capaz de leer este arbol y ejecutar un ciclo sobre cada evento contenido en el arbol. Ahora nos preguntamos: ¿qué sucede si me dan un arbol cualquiera con datos y necesito hacer un análisis? ¿Como proceder?

En realidad esta pregunta ya la hemos contestado en la Nota Suplementaria N1 del TC2. Podemos hacer que ROOT cree una nueva clase basada en cualquier arbol que le pasemos y sobre esta clase construir nuestro análisis.

Este ejercicio se basa en un arbol que contiene un numero pequeño de eventos de datos reales obtenidos por el experimento CMS del LHC, durante la corrida del 2011. Para aquel que trabaja en CMS y desee conocer como se extrajeron estos datos, ver la siguiente Nota suplementaria: EPUniandes2012Ntupler.

• Hagamos brevemente una exploracion del archivo real_data_sample1.root, el cual debe encontrarse en su carpeta, usando el TBrowser de ROOT:

root -l real_data_sample1.root
root [0] new TBrowser

• Vemos el arbol con los datos que queremos analizar. Tomar nota del nombre tal como aparece: t1.

Show results... Hide results...

• Este arbol se construyo exclusivamente con informacion parcial de muones que han sido reconstruidos por el detector CMS. En realidad se escogieron eventos con 2 o mas muones por evento. Se guardaron las siguientes variables con propiedades cinematicas de estos:
  • muons_fPx: componente en X del momento
  • muons_fPy: componente en Y del momento
  • muons_fPz: componente en Z del momento
  • muons_fE: Energia
  • muons_fC: carga del muon
  • muons_fEta: pseudorapidez
  • nMuons: numero de muones en el evento

• Permaneciendo en ROOT, ahora procedemos de forma similar a lo sugerido en la nota Suplementaria N1 del TC2:

root [1] _file0->cd();
root [2] TTree *tree = (TTree*)gDirectory->FindObjectAny("t1");
root [3] tree->MakeClass("DataAnalysis");
Info in : Files: DataAnalysis.h and DataAnalysis.C generated from TTree: t1
root [4] .q

• Salir de ROOT y abrir los dos archivos creados:

gedit DataAnalysis.h DataAnalysis.C

• Explorar: el primer archivo, esto es el .h contiene la declaración de la clase DataAnalysis y la definición de algunos métodos que ROOT implementa para la lectura del árbol de datos. El segundo archivo, es decir el .C contiene la definición del método Loop(). Este es el método en cargado de hacer un ciclo sobre las entradas del árbol. Mantengan abiertos los archivos, los vamos a editar en la siguiente parte.

Ejercicio 2 - Preparativos

En este ejercicio, vamos a implementar algunos cambios a la clase DataAnalysis. El primer cambio será adicionar un constructor a la clase. Para ello editar el archivo DataAnalysis.h (si no está abierto, abrirlo):

gedit DataAnalysis.h

• Necesitaremos incluir las declaraciones de la clase que posee ROOT para la manipulación de 4-vectores. Por ello, justo al final del ultimo "#include", tendremos que adicionar lo siguiente:

#include < TLorentzVector.h >

• Ahora, necesitamos bajar en el código hasta aproximadamente la linea No 40 y copiar la siguiente linea:

DataAnalysis(TTree *tree=0); // este es el constructor normal - no cambiar - insertar las lineas de abajo
// Nuevo constructor - toma como argumento el nombre del archivo que vamos a analizar
DataAnalysis(const char *); 

Tal como lo dice el comentario, esta linea nos permite crear un objeto de tipo DataAnalysis pasando como argumento el archivo a analizar. Esto es muy conveniente. El constructor base que ROOT le da a esta clase, toma como argumento un TTree y por consiguiente tendríamos que abrir por separado nuestro archivo de datos, extraer el árbol y pasarlo a DataAnalysis.

• Guardar. Ahora tenemos que editar el archivo DataAnalysis.C y terminar de definir este constructor (si no está abierto, abrirlo).

gedit DataAnalysis.C

#define DataAnalysis_cxx
#include "DataAnalysis.h"
#include < TH2.h >
#include < TStyle.h >
#include < TCanvas.h >

//Insertar las siguientes lineas  - esta es la definicion del nuevo constructor

 DataAnalysis::DataAnalysis(const char * fname)
{

  // Cambiar el constructor - preferible que tome como argumento un nombre de archivo
  TTree * tree; //Declaramos un apuntador de tipo TTree
  TFile * f = TFile::Open(fname); //aqui abrimos el archivo con nombre fname
  f->cd();
  tree = (TTree*)gDirectory->Get("t1"); //extraemos un apuntado al arbol con nombre t1
  Init(tree); //Incializamos la conexion en memoria con este t1
  
} 

• Perfecto! Ahora estamos listos para proceder tal como hicimos en los Talleres anteriores, en donde modificamos el metodo Loop() y analizar los datos.

Ejercicio 3 - Análisis de datos

• En este ejercicio haremos un histograma de masa invariante de un par muones de cargas opuestas. La información necesaria se encuentra contenida en el archivo de datos. Nuevamente, trabajamos en el método Loop() que se encuentra en el archivo DataAnalysis.C. Dado que esto es código que ROOT ha creado automáticamente, aparecen toda clase de comentarios que ustedes pueden leer (describen como funciona esta clase) y luego borrar o ignorar.

• En este ejercicio, el interés es hacer un histograma de masa invariante, por lo tanto lo primero que debemos hacer es definir el histograma. Podemos hacerlo en el rango siguiente XMIN=0 XMAX = 120.0 y un numero de bins de 120. Intente hacerlo sin mirar la respuesta -recuerde que debe ir justo antes del inicio del ciclo sobre eventos:

Show results... Hide results...

 //Mi primer histograma de masa invariante con datos REALES del LHC!!!
 TH1F * h_DiMuonInv = new TH1F("mumuInv","Masa invariante de muones", 120, 0, 120); 

• Dado que ahora trabajaremos con datos, vamos a darle un tratamiento diferente a cada punto en histograma, por lo que adicionaremos la opción siguiente para que ROOT trate correctamente los errores:

//Mi primer histograma de masa invariante con datos REALES del LHC!!!
TH1F * h_DiMuonInv = new TH1F("mumuInv","Masa invariante de muones", 120, 0, 120);
h_DiMuonInv->Sumw2();

• Lo que sigue debe ser muy familiar, pues lo hemos hecho en los Talleres pasados. Un dificultad es que los eventos poseen 2 o mas muones y por lo tanto necesitamos formar todas las parejas posibles. Ademas, vamos a una necesitar una condición que seleccione parejas de muones de cargas opuestas. Mi propuesta es la siguiente (usted puede pensar en una alternativa):

 // Estamos dentro del ciclo sobre eventos - despues de if (Cut(ientry) < 0) continue

  //Hacer un loop sobre los muones - necesitamos definir dos 4-vectores para el muon 1 y el muon 2

  TLorentzVector muonOne;
  TLorentzVector muonTwo;

  // El maximo numero de muones en este evento, lo tenemos gracias a la variable nMuons, contenida en el TTree
  int maxMuons = nMuons;

  // Iniciamos un ciclo con indice i
    
    for( int i = 0; i < maxMuons; ++i) {

      // Necesitamos formar parejas de muones i-j
      //   luego puedo iniciar un ciclo secundario sobre el indice j y que empieze en i+1

      for( int j = (i+1); j < maxMuons; ++j) {

        // condicion: si el indice i es diferente del indice j  -Y-  la carga de los muones es opuesta entonces proceder 
        if ( ( i != j ) && (muons_fC[i] != muons_fC[j] )) {
                    
        } // cerrar la condicion
                
      } // cerrar ciclo sobre las j
      
    } // cerrar cicclo sobre las i
   

• Bien, tenemos ahora el esqueleto de nuestro análisis basado solo en los muones reconstruidos. A diferencia de las simulaciones aqui no manejamos codigos. Confiamos en que el detector CMS, su sistema de adquisición de datos y de reconstrucción nos entreguen los mejores candidatos a muones (en realidad asi es, la eficiencia esta por encima de 95%). Ahora solo falta adicionar dentro de la condición en la que se forman las parejas de muones de carga opuesta, los 4-vectores de los dos muones, combinarlos, calcular su masa invariante y llenar el histograma.

 // Esto va dentro de la condicion
          muonOne.SetPxPyPzE( muons_fPx[i], muons_fPy[i], muons_fPz[i], muons_fE[i] ); // muon 1 va con el indice i
          muonTwo.SetPxPyPzE( muons_fPx[j], muons_fPy[j], muons_fPz[j], muons_fE[j] ); // muon 2 va con el indice j

          TLorentzVector combVector = (muonOne + muonTwo);
          float invMass = combVector.M();  
          h_DiMuonInv ->Fill(invMass); 

• En las dos primeras lineas, establecemos los 4-vectores para muon 1 y muon 2, luego los sumamos y calculamos la masa invariante del 4-vector resultante, gracias al metodo M(). Finalmente llenamos el histograma.

• Tal como hicimos en ejercicios anteriores, despues de finalizar el ciclo sobre eventos, damos la instruccion de graficar la distribucion de masa invariante de los dos muones.

 //graficar histogramas
h_DiMuonInv ->Draw(); 

• Guardar y cerrar.

• Ahora para correr este codigo de forma similar a como hicimos en el TC2:

root -l
root [0] gROOT->LoadMacro("DataAnalysis.C")
(Int_t)0

• Si todo sale bien, no deberia haber ningun mensaje de error. Si aparece alguno, hay que leer el mensaje. Posiblemente ha quedado algo mal escrito durante la modificacion realizadas (revisar).

• Ahora para correr este codigo, necesitamos dos lineas. La primera sera crear un apuntador a un objeto de la clase DataAnalysis y luego a este apuntador, hacer que corra el metodo Loop(). El constructor de un objeto DataAnalysis toma como argumento el nombre del archivo (va entre comillas) en donde se encuentran los datos simulados.

root [1] DataAnalysis * data = new DataAnalysis("real_data_sample1.root");
root [2] data->Loop()

• Voila! El resultado es fantastico: podemos ver claramente un pico alrededor de los 90 GeV/c^2. ¿Sabe usted a que corresponde?

Show results... Hide results...

En resumen: Si recibimos un archivo con datos que se encuentren con una estructura arbol en ROOT, sabemos como atacar el problema de analizarlos utilizando las herramientas que ROOT no provee. Ahora bien, este no es el fin de la historia y hacer un analisis sobre datos finales tiene muchas complejidades.

• La colaboracion CMS posee una grafica oficial con L=1.1 fb-1 de datos del 2011 (muy similar a la cantidad de datos en los que basamos este sencillo analisis).

Show results... Hide results...

• Si usted quisiera ver si en estos eventos se pueden reconstruir J/Psi's, como modificaría el código?

Eso seria todo. Si tienen inquietudes, preguntas o dudas pueden plantearlas.

-- AndresOsorio - 27-May-2012