ANÁLISIS DE DATOS PÚBLICOS DE CMS CON R

Primero hay que cargar un archivo con los datos que queramos utilizar

In [2]:
dimuon  <- read.csv("https://indico.cern.ch/event/577111/contributions/2659823/attachments/1498147/2337502/MuRun2010B.csv")

Y ahora vamos a hacer un histograma con las masas invariantes de los dimuones de signos opuestos

In [3]:
hist(dimuon$M, breaks = 108, main = "Espectro de masas invariantes de dimuones de signos opuestos", xlab = "M / GeV", ylab = "Número de sucesos por bin")

Los bines del histograma tienen 1 GeV, es decir, las masas están agrupadas en clases de 1 GeV de ancho. De momento los elegimos mediante prueba y error...

Parece que hay algo de estructura entre 0 y 10 GeV y, quizá, habría que mirar cerca de 90 GeV. Desde luego, esa zona a la derecha de 30 GeV hay que verla con otra escala.

Vamos a mirar primero el histograma sólo entre 0 y 10 GeV

In [4]:
hist(dimuon$M, breaks = 108, main = "Espectro de masas invariantes de dimuones de signos opuestos", xlab = "M / GeV", ylab = "Número de sucesos por bin", xlim = c(2, 15))

Claramente los bines son demasiado anchos, pasemos (¡prueba y error!) a bines 100 veces más pequeños, de 10 MeV

In [5]:
hist(dimuon$M, breaks = 10800, main = "Espectro de masas invariantes de dimuones de signos opuestos", xlab = "M / GeV", ylab = "Número de sucesos por bin", xlim = c(2, 110))

Ahora vamos a ampliar ese pico en torno a los 3 GeV

In [6]:
hist(dimuon$M, breaks = 10800, main = "Espectro de masas invariantes de dimuones de signos opuestos", xlab = "M / GeV", ylab = "Número de sucesos por bin", xlim = c(2, 12))

Pues ahí, con una masa de 3.1 GeV hay un pico que se puede interpretar como que hemos redescubierto una partícula, la J/psi. También hay otras regiones que merece la pena explorar en detalle (digamos que entre 3,5 y 4 GeV y también entre 9 y 10 GeV...)

Ahora vamos a modificar el histograma pasando de barras a puntos (luego habría que poner barras de error. etc).

In [7]:
# Primero extraemos los datos del histograma. Ahora bines de 20 MeV
datos <- hist(dimuon$M, breaks = 5400, plot = F)
In [8]:
# Y ahora los representamos: datos$mids son los centros de los bines y datos$counts el número de eventos de cada bin
plot(datos$mids, datos$counts, pch = 16, main = "Espectro de masas invariantes de dimuones de signos opuestos", xlab = "M / GeV", ylab = "Número de sucesos por bin")

Vamos a ampliar un poco la región de masas entre 9 y 10 GeV (bueno, GeV/c^2)

In [13]:
plot(datos$mids, datos$counts, pch = 16, main = "Espectro de masas invariantes de dimuones de signos opuestos", xlab = "M / GeV", ylab = "Número de sucesos por bin", xlim = c(9, 10), ylim = c(50, 275))

El siguiente paso lógico sería poner barras de error, arreglar los ejes... pero lo dejaremos para otra ocasión.

Y para concluir, volmemos al rango completo y pasamos a una doble escala logarítmica

In [17]:
plot(log10(datos$mids), log10(datos$counts), pch = 20, cex = 0.4, main = "Espectro de masas invariantes de dimuones de signos opuestos", xlab = "log M / GeV", ylab = "log N sucesos por bin")

Ahora lo suyo sería poner barras de error y rotular bien los ejes...