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:


tar xzvf ejercicios-Parte3.tgz
cd ejercicios-Parte3

En esta carpeta, se encuentra el codigo necesario para la hacer los ejercicios. Tambien existe una carpeta solucion que contiene el codigo completo.

Ejercicio 1 – Creación de un simple Histograma

El primer ejercicio consistirá en la creación de un Histograma con informacion de la cinematica de algunas particulas. Editemos el script que se trabajo al final del Taller 2 y que simula eventos de Supersimetria:

gedit pythia8_susy.C

  • Quitar el comentario en la linea siguiente del codigo y escribir el numero de bins 50 y el rango del histograma XMIN=0 y XMAX=100:

//TH1F * h1 = new TH1F("hpt", "Muons Pt", 50, 0, 100);

  • Esta linea define un Histograma 1D ( clase TH1F ).

  • Bajar hasta donde está el ciclo de generación y buscar la linea con el comentario "// Particle Loop" y agregrar el siguiente codigo:

while( iev < maxEvts ) {
    
    pythia8->GenerateEvent();
    if (iev < 1) pythia8->EventListing();
    pythia8->ImportParticles(particles,"All");
    
    Int_t np = particles->GetEntriesFast();

    // Particle loop
    
    for (Int_t ip = 0; ip < np; ip++) {
      
      //Leemos informacion de cada particula

      TParticle* part = (TParticle*) particles->At(ip); //Esta es la particula i
      Int_t ist = part->GetStatusCode(); //obtenemos su estatus
      Int_t pdg = part->GetPdgCode(); // y su codigo PDG
      
       //1. Introducir una codicion para separar los muones

      //Queremos graficar el momento Transverso de los muones Pt == Sqrt(px*px + py*py)
      //2. Obtener aqui las componentes Px y Py del momento de las particulas

      //3. Llenar el histograma h1 con el momento transversar

      
    } // No olvidar! cerrar el loop que hemos introducido 

Estudiemos que esta pasando aqui:

  • El primer ciclo, es un ciclo sobre la generacion de eventos . Por ello es que la condicion es hacer este ciclo hasta que alcance el maximo numero de eventos que hemos especificado desde el principio del script. La generacion esta instruida por el comando pythia8->GenerateEvent();

  • El segundo ciclo, aquel que hemos introducido, hace por cada evento un ciclo sobre el numero maximo de particulas que Pythia ha generado. El maximo numero de particulas se obtiene con el metodo Int_t np = particles->GetEntriesFast();

  • Los muones tienen por código de identificación 13 y -13 (para el anti-mu). Volvemos a este comentario mas adelante

  • Cada partícula está representada por un objeto de tipo TParticle y queremos graficar el momento transversal de cirtas particulas
    • Mirando en la documentación TParticle, como podríamos obtenerlo?

  • Por el momento, el codigo introducido no hace aun lo que queremos. Recordemos, queremos hacer 3 cosas:

  1. Entre la particulas, identificar aquellas que son muones -es decir con codigo 13 y -13
  2. Un vez identificadas, tenemos calcular el momento transversal
  3. Llenar el histograma que creamos para hacer

  • Por lo tanto, debemos copiar el siguiente codigo - guiarse por los comentarios:


      //1. Introducir una codicion para separar los muones
      if ( pdg == 13 || pdg == -13 ) {

      //Queremos graficar el momento Transverso de los muones Pt == Sqrt(px*px + py*py)
      //2. Obtener aqui las componentes Px y Py del momento de las particulas

      float px = part->Px();
      float py = part->Px();
      float Pt = sqrt(px*px + py*py);

      //3. Llenar el histograma h1 con el momento transversar
      h1->Fill( Pt );
 
      } ///no olvidar! cerrar el condicional 

* Bien! Esto nos sirve para llenar el histograma. Ahora, si queremos visualizar el histograma debemos introducir al final de todo los ciclos lo siguiente:

    // Ahora podemos dibujar el histograma y que lo muestre en pantalla
    h1->Draw();

  • Guarde las modificaciones y cierre. Ejecutar ahora el script desde la terminal:

root -l pythia8_susy.C

La opción -l ("menos L minuscula") evita que salga la imagen de saludo de ROOT, lo que nos ahorra algunos milisegundos.

muonpt-1.jpg

Ejercicio 2 – PythiaAnalsis – Una clase para hacer el análisis sobre datos simulados

Ahora lo que queremos es introducir una técnica para el análisis de datos simulados.

Hemos trabajadoe en el script pythia8_susy.C, el cual genera los eventos y los guarda en un archivo de formato de ROOT (manejado por la clase TFile). En el script, la creacion del archivo ocurre en la siguiente linea:

//Definir archivo de salida
TFile * outfile = new TFile("eventos_pythia8_SUSY.root","RECREATE");

Dentro de este archivo existe lo que se conoce como un TTree – un árbol – en el que cada rama corresponde a las variables que poseen las partículas y las hojas son los valores que toman por evento. Este árbol se define en la siguiente linea del código:

TTree*tree= new TTree("tree","Arbol con particulas segun Pythia8");
tree->Branch("particles",&particles);

¿Cómo hacer un analisis de los datos contenidos en este árbol? Lo haremos utilizando una clase llamada PythiaAnalysis. Esta clase es automáticamente generada por ROOT basada en la estructura que tiene árbol (ver nota Suplementaria No 1). La clase PythiaAnalysis se declara y define en dos archivos:

  • PythiaAnalysis.h: declaración de la clase, sus variables y metodos
  • PythiaAnalysis.C: definición del metodo Loop. Este metodo realiza un ciclo/loop sobre los eventos contenidos en el árbol

Miremos el archivo PythiaAnalysis.h

gedit PythiaAnalysis.h

  • El este archivo se declara la clase PythiaAnalysis. Contiene toda una serie de arreglos, los cuales definen el contenido por evento de cada particula simulada.
  • Tambien contiene unos metodos (en el contexto de C/C++). Entre los mas importantes, se encuentran un constructor y otros metodos que se encargan de asociar y llenar cada uno de de los arreglos con los datos contenidos en el arbol.
  • Para nosotros, el metodo que vamos a modificar es el llamado Loop, el cual hace un ciclo sobre los eventos generados. Este metodo se define en el archivo PythiaAnalisis.C. Editemos entonces este archivo:

gedit PythiaAnalysis.C

Intentemos hacer de nuevo un histograma que nos cuente el numero de muones por evento. Nos concentraremos en implementar el método Loop( ) que se encarga de hacer la iteración sobre los eventos almacenados.

  • Para ello debemos adicionar un segundo ciclo, que nos itere sobre la particula i-esima, parecido a lo que hicimos en el ejercicio anterior:
  • En este ciclo, preguntamos por el PDG_ID de la particula y comparamos
  • Adicionamos una condición y un contador para muones
  • Finalizado este loop podemos entonces graficar la distribución del número de muones por evento

Nota: la lectura del TTree nos hace perder los objetos partículas como tales (TParticle). Su información se encuentra ahora en forma de arreglos o arrays a la C/C++ (ver PythiaAnalysis.h).

  • Uno tiene dos opciones:

  1. leer y manipular los arreglos
  2. leer los arreglos y crear de ellos objetos TParticle

Para este ejercicio, haremos una sencilla modificacion del metodo Loop() y utilizaremos los arreglos.


   //Adicionar un histograma para contar el numero de muones
  TH1F * h1 = new TH1F("nmu","Muones por eventos", 50, 0, 15); 

  // Empieza loop sobre eventos  
  for (Long64_t jentry=0; jentryGetEntry(jentry);   nbytes += nb;
    
    //Implementar aqui el analisis
    
    int maxpart = particles_; //maxpart contiene el numero maximo de particulas
    int n_muons = 0; //contador para el numero de muones _por_ evento
    
    for (int i = 0; i < maxpart; ++i ) // ciclo sobre particulas por evento
    {
      
      if( abs( particles_fPdgCode[i] ) == 13 ) //condicion que nos identifica los muones ( de codigo +13/-13)
        ++n_muons; //incrementamos el contador
      
    }
    
    h1->Fill(n_muons); //llenamos el histograma 


  }//cierra loop sobre eventos  
   
  //Graficar el histograma del numero de muones por evento
  h1->Draw(); 

Como correr este codigo? Un forma sencilla es la de cargar esta clase al inicio de la sesion en ROOT:

root -l
root [0] gROOT->LoadMacro("PythiaAnalysis.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 realizada al metodo Loop().

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

root [1] PythiaAnalysis * susy = new PythiaAnalysis("eventos_pythia8_SUSY.root");
root [2] susy->Loop()

  • En este ejemplo, susy es un apuntador a un objeto de tipo PythiaAnalysis. Una vez creado el objeto al que apunta, podemos llamar el metodo Loop y hacer el analisis sobre los eventos simulador.

nmuones-2.png

Para hacer esto de una forma mas automatica, ver nota Suplementaria No 2.

Ejercicio 3 – JetPythiaAnalysis – Una clase para hacer el análisis sobre datos simulados incluyendo Jets

Ahora queremos introducir una técnica para el análisis de datos simulados usando los algoritmos de reconstrucción ded Jets que nos da el paquete FastJet (incluido en el Live-CD que ustedes tienen). Esto lo hemos implementado en una clase llamada JetPythiaAnalysis la cual posee las mismas caracteristicas de PythiaAnalysis, pero además incluye la posibilidad de reconstruir Jets.

  • Abrir los archivos JetPythiaAnalysis.h y JetPythiaAnalysis.C : explorar
    • Al igual que en el caso anterior, en JetPythiaAnalysis.h se declara la clase (con sus variables y metodos)
    • JetPythiaAnalysis.C contiene la definicion del metodo Loop()
  • Vamos a editar el archivo JetPythiaAnalysis.C

La idea es la de construir Jets a partir de partículas estables o cuyo tiempo de vida es suficiente para que alcancen a llegar a los distintos instrumentos de deteccion que se colocan en un detector de particulas. Por ejemplo, dentro de este tipo de particulas tendriamos muones, electrones, piones, kaones. En este ejercicio haremos la reconstruccion de Jets segun un algorithmo en especial, el Kt-jet y graficaremos la distribucion de momento transversal. Para ello los pasos a seguir son:

  1. Definir el algoritmo para la reconstruccion de Jets
  2. Definir el histograma
  3. Hacer un ciclo sobre particulas por evento y seleccionar aquellas que sean estables y que detectables. Por detectables nos referimos a que particulas que podrian ser vistas en el detector y no particulas que como neutrinos escaparian a la deteccion.
  4. Llenar un contenedor con estas particulas
  5. Correr el algorithmo de reconstruccion de Jets
  6. Hacer un ciclo sobre los Jets resultantes, extraer su momento transversal y pasar el valor al histograma
  7. Una vez terminado el ciclo sobre eventos, graficar

Editar el archivo JetPythiaAnalysis.C siguiendo los comentarios que alli se encuentran:

gedit JetPythiaAnalysis.C

e introducir las siguientes lineas de codigo:

  • Definir el algoritmo para reconstruir jets

  //Definir el tipo de Algoritmo para reconstruir Jets y pasarle el parametro requerido
   double R = 0.7;
   fastjet::JetDefinition jet_def(fastjet::kt_algorithm, R);

  • Definir el histograma que vamos a llenar

  //Definir un histograma para el momento transversal de los Jets
 TH1F * h_JetPt = new TH1F("JetsPt","Momento transversal de Jets ",50, 0, 100);

  • En el ciclo de eventos, insertar un ciclo sobre particulas y colocar alli las condiciones necesarias para seleccionar las particulas estables y detectables y llenar un contenedor.

    int np = 0;
     int max_part = particles_; //Maximo numero de particulas

    // Definit un contenedor de Particulas para pasarle a FastJet
     std::vector particles;
    
    //Implementar aqui el analisis
     while ( np < max_part )  {
      
      //llenar aqui con las particulas estables: estas particulas serian aquellas que no tienen hijas

      if ( (particles_fDaughter[np][0] == -1) && (particles_fDaughter[np][1] == -1) ) 
      {

        bool isDetectable = true;

        // tenemos que saltarnos aquellas particulas que aun siendo estables, no serian detectables directables
	// 12 = nu_e ; 14 = nu_mu ; 16 = nu_tau ; 10000022 = SUSY LSP / neutrinalino

        if ( abs (particles_fPdgCode[np]) == 12 || 
	     abs (particles_fPdgCode[np]) == 14 || 
	     abs (particles_fPdgCode[np]) == 16 || 
	     abs (particles_fPdgCode[np]) == 10000022 )
			isDetectable = false;

        // lenamos el contenedor con particulas estables
        if( isDetectable ) particles.push_back( fastjet::PseudoJet( particles_fPx[np], 
                          	                                    particles_fPy[np],
                          	                                    particles_fPz[np],
                          	                                    particles_fE[np] ) );
      
       } //Este cierra la condicion
      
      ++np; // pasamos a la siguiente particula
       
    }

  • Ejecutar el algoritmo de reconstruccion de Jets y hacer un ciclo sobre el resultado. Extraer el momento transversal de cada Jet y llenar el histograma con este valor.

    if ( particles.size() <= 0 )
      continue;
    
    // Corremos ahora FastJet: hacer sobre las particulas la identificacion de Jets
     fastjet::ClusterSequence Cluster(particles, jet_def);

    // Objeto Cluster contiene los jets: podemos guardarlos en un contenedor de Jets
     std::vector jets = Cluster.inclusive_jets();
    
    // Ahora podemos hacer un ciclo sobre los jets reconstruidos y extraer el momento transversal
    // - llenar histograma

     for (unsigned i = 0; i < jets.size(); i++) {

      h_JetPt->Fill( jets[i].perp() );
      
    }

     jets.clear(); //Limpiar el contenedor de jets y particulas en preparacion para el proximo evento
     particles.clear(); 

  • Por fuera del ciclo sobre los eventos, graficar el histograma:

  // Dibujar la distribucion de momento transversal
   h_JetPt->Draw();

  • Listo, esas serian todas las modificaciones a hacer al codigo. Para correr vamos a necesitar algunas cosas mas. Editar el archivo siguiente y seguir comentario:

gedit rootlogon.C

   //Para el ejercicio 3: needCompilation = true
   bool needCompilation =  true;

  • Que ocurre? Resulta que para poder usar FastJet dentro de ROOT, necesitamos una opcion que cargue la libreria de FastJet y que compile el codigo que hemos realizado, a diferencia de los pasos anteriores. Esto lo podemos instruir para que se haga automaticamente al inicio de ROOT desde este archivo rootlogon.C (ver nota Suplementaria N2).

  • Si todo ha quedado bien escrito, ROOT compila exitosamente las dos clases PythiaAnalysis y JetPythiaAnalysis.

root -l
Info in : creating shared library /ejercicios-Parte3/solucion/./PythiaAnalysis_C.so
Info in : creating shared library /ejercicios-Parte3/solucion/./JetPythiaAnalysis_C.so
root [0]

  • Para correr JetPythiaAnalysis, tenemos que hacer dos pasos muy similares a los del ejercicio N2:

root -l
Info in : creating shared library /ejercicios-Parte3/solucion/./PythiaAnalysis_C.so
Info in : creating shared library /ejercicios-Parte3/solucion/./JetPythiaAnalysis_C.so
root [0] JetPythiaAnalysis * jets = new JetPythiaAnalysis("eventos_pythia8_SUSY.root");
root [1] jets->Loop(); 
#--------------------------------------------------------------------------
#                      FastJet release 2.4.2
#            Written by M. Cacciari, G.P. Salam and G. Soyez            
#                         http://www.fastjet.fr                         
#								      	   
# Longitudinally invariant Kt, anti-Kt, and inclusive Cambridge/Aachen  
# clustering using fast geometric algorithms, with area measures and optional
# external jet-finder plugins.                                          
# Please cite Phys. Lett. B641 (2006) [hep-ph/0512210] if you use this code.
#								      	   
# This package uses T.Chan's closest pair algorithm, Proc.13th ACM-SIAM 
# Symp. Discr. Alg, p.472 (2002), S.Fortune's Voronoi algorithm and code .
#-------------------------------------------------------------------------

  • El mensaje en pantalla es un mensaje de bienvenida de FastJet. El resultado de nuestro analisis sera algo como esto:

jets-Pt-3.png

Eso sería todo! Buena suerte.

Nota Suplementaria N1

Estos son los pasos para crear una clase para hacer el analisis de un TTree:

root -l
root [0] TFile *_file0 = TFile::Open("eventos_pythia8_SUSY.root");
root [1] _file0->cd();
root [2] TTree *t1 = (TTree*)gDirectory->FindObjectAny("tree");
root [3] t1->MakeClass("MyFirstROOTclass");
Info in : Files: MyFirstROOTclass.h and MyFirstROOTclass.C generated from TTree: tree

Nota Suplementaria N2

Cada vez que uno inicia ROOT, el programa busca en la carpeta en donde se encuentra parado un script llamado rootlogon.C. Si no existe, ROOT inicia con los valores que trae por por defecto. Si existe, entonces ejecuta este script, dando la oportunidad de que el usuario introduzca cambios en las configuracion inicial. Nosotros podemos usar este script para que carge automaticamente las clases que necesitamos. Este seria un ejemplo del contenido de dicho script:

// Contenido del archivo rootlogon.C
void rootlogon() 
{
gROOT->LoadMacro("PythiaAnalysis.C");
}

De esta forma nos evitamos este paso. Para ejecutar, el Loop sobre unos datos podemos guardar las dos lineas en un macro. Un macro almacena los comandos que dariamos a ROOT en una sesion interactiva, los cuales deben ir entre dos corchetes { y }. Al igual que en C/C++ cada comando debe ir separado por un ;. Por ejemplo, uno podria almacenar todo en el siguiente macro:

// Contenido del archivo runAnalysis.C
{

PythiaAnalysis * susy = new PythiaAnalysis("eventos_pythia8_SUSY.root");
susy->Loop();

}

Este macro se puede ejecutar de dos formas:

  • Desde linea de comando:

root -l runAnalysis.C

  • Desde ROOT con el comando .x (punto x):

root -l 
root [0] .x runAnalysis.C

-- AndresOsorio - 15-May-2012

Topic attachments
I Attachment HistorySorted ascending Action Size Date Who Comment
Compressed Zip archivetgz ejercicios-Parte3.tgz r1 manage 49.9 K 2012-05-24 - 18:19 AndresOsorio material P3
PNGpng jets-Pt-3.png r1 manage 9.0 K 2012-05-24 - 17:32 AndresOsorio Jetspt-3
JPEGjpg muonpt-1.jpg r1 manage 39.3 K 2012-05-22 - 17:10 AndresOsorio muonpt-1
PNGpng nmuones-2.png r1 manage 8.5 K 2012-05-22 - 18:46 AndresOsorio nmuones-2
Edit | Attach | Watch | Print version | History: r11 < r10 < r9 < r8 < r7 | Backlinks | Raw View | WYSIWYG | More topic actions
Topic revision: r11 - 2013-08-22 - AndresOsorio
 
    • 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