Skip to content
Snippets Groups Projects
Commit 5a6aed1e authored by Luis Pascual's avatar Luis Pascual
Browse files

Upload New File

parent cce278d4
Branches master
No related tags found
No related merge requests found
%% Cell type:markdown id:racial-confidence tags:
# Ajuste de los parámetros de un modelo a datos "toy"
### Este notebook está organizado de la siguiente forma:
1. Definición de modelos de señal y fondo para generar nuestros "toys"
2. Visualización de los datos
3. Ajuste de los parámetros minimizando el $\chi^{2}$ entre el modelo y los datos.
4. Profundizando en el proceso de minimización
%% Cell type:markdown id:decreased-television tags:
Antes de empezar, importamos las librerías necesarias
%% Cell type:code id:equivalent-toyota tags:
``` python
import ROOT
import numpy as np
```
%% Output
Welcome to JupyROOT 6.24/02
%% Cell type:markdown id:handy-order tags:
Creamos lienzos sobre los que mostrar los plots
%% Cell type:code id:committed-bench tags:
``` python
cModels = ROOT.TCanvas('cModels','cModels',1200, 400) #este lo usaremos para los modelos
cModels.Divide(3,1) #lo dividimos en una red the 3x1
cToys = ROOT.TCanvas('cToys','cToys',1200, 600) #este lo usaremos para los modelos
cToys.Divide(2,1) #lo dividimos en una red the 3x1
c1 = ROOT.TCanvas('c1','c1',900, 600)#este es comodin, por si queremos dibujar algo
ROOT.gStyle.SetPalette(1) #opcional: cambia la paleta de colores en los plots.
```
%% Cell type:markdown id:genetic-incidence tags:
## Modelos de señal y fondo y "toys"
El modelo que vamos a construir está "inspirado" en el descubrimiento del boson de Higgs en su canal de desintegración difotón:
- La forma de la señal es descrita por una distribución Gaussiana.
- La forma del fondo es descrita por una distribución exponencial.
El número de eventos totales se podría expresar de la siguiente forma:
$ N_{total} = N_{sig}\times PDF_{sig} + N_{bkg}\times PDF_{bkg} = N_{total} \left(f_{sig}\times PDF_{sig} +(1-f_{sig})\times PDF_{bkg}\right)$
donde:
- $N_{total,sig,bkg}$ son el número de eventos totales, de señal y de fondo respectivamente.
- $PDF_{sig}$ y $PDF_{bkg}$ son las funciones de densidad de probabilidad de la señal y fondo respectivamente.
- $f_{sig}$ es la fracción de eventos de señal con respecto al número de eventos totales $\left(=N_{sig}/(N_{sig}+N_{bkg})=N_{sig}/N_{total}\right)$
%% Cell type:markdown id:clean-sunset tags:
Primero, definimos las variables y parámetros que necesitamos:
%% Cell type:code id:yellow-nation tags:
``` python
fSig = 0.05
nEntries = 10000
minVal = 100.000
maxVal = 160.000
```
%% Cell type:markdown id:shared-journey tags:
Cada componente del modelo (señal o fondo) va a ser definido como una PDF en el intervalo definido arriba (minVal,maxVal).
%% Cell type:markdown id:academic-prospect tags:
Como hemos dicho, la señal está descrita por una distribución gaussiana, centrada en un valor de masa 'mass' y con una anchura 'sigma'.
En este caso tenemos herramientas que normalizan automáticamente la distribución.
%% Cell type:code id:lesbian-ordinance tags:
``` python
#PDF de la señal, más info: https://root.cern.ch/root/html524/TMath.html#TMath:Gaus
signalModel = ROOT.TF1('signalModel','TMath::Gaus(x,[0],[1],1)',minVal,maxVal)
signalModel.SetParNames('mass','sigma') #nombres de los parámetros
signalModel.SetParameters(125,2.4) #valores de los parámetros
cModels.cd(1)
signalModel.Draw()
```
%% Cell type:markdown id:cleared-scottish tags:
El fondo, descrito por una exponencial, require una normalización explícita.
Nosotros hoy la vamos a hacer utilizando directamente la integral de la función calculada de forma numérica.
Sin embargo, en este caso tiene solución análitica (es un buen ejercicio, comprobad que la sabéis obtener...)
%% Cell type:code id:hollow-dispute tags:
``` python
#PDF del fondo, una exponencial
bkgModel = ROOT.TF1('bkgModel','TMath::Exp(-(x-[1])/[0])',minVal,maxVal)
bkgModel.SetParameters(50,100)
normBkg = 1./bkgModel.Integral(minVal,maxVal)
print('Integral fondo en el rango utilizado: {}'.format(bkgModel.Integral(minVal,maxVal)))
bkgModel = ROOT.TF1('bkgModel','{}*TMath::Exp(-(x-[1])/[0])'.format(normBkg),minVal,maxVal)
bkgModel.SetParNames('tau','delta')
bkgModel.SetParameters(50,100)
cModels.cd(2)
bkgModel.Draw()
```
%% Output
Integral fondo en el rango utilizado: 34.940289404389894
%% Cell type:markdown id:connected-march tags:
Ahora, combinamos ambas PDF en una sola, incluyendo la cantidad de eventos de señal y fondo.
%% Cell type:code id:composed-hunter tags:
``` python
fullModel = ROOT.TF1('fullModel','(signalModel*[0]+[1]*bkgModel)',minVal,maxVal)
fullModel.SetParName(0,'nsig')
fullModel.SetParName(1,'nbkg')
fullModel.SetParameter('nsig',nEntries*fSig)
fullModel.SetParameter('nbkg',nEntries*(1-fSig))
fullModel.FixParameter(2,100)
fullModel.Print('v')
cModels.cd(3)
fullModel.Draw()
```
%% Output
Formula based function: fullModel
fullModel : (signalModel*[0]+[1]*bkgModel) Ndim= 1, Npar= 6, Number= 0
Formula expression:
((TMath::Gaus(x,[mass],[sigma],1))*[nsig]+[nbkg]*(0.028620255213866665*TMath::Exp(-(x-[delta])/[tau])))
List of Variables:
Var 0 x = 0.000000
List of Parameters:
Par 0 nsig = 500.000000
Par 1 nbkg = 9500.000000
Par 2 delta = 100.000000
Par 3 mass = 125.000000
Par 4 sigma = 2.400000
Par 5 tau = 50.000000
Expression passed to Cling:
#pragma cling optimize(2)
Double_t TFormula____id14889634374136314836(Double_t *x,Double_t *p){ return ((TMath::Gaus(x[0],p[3],p[4],1))*p[0]+p[1]*(0.028620255213866665*TMath::Exp(-(x[0]-p[2])/p[5]))) ; }
%% Cell type:code id:paperback-recommendation tags:
``` python
cModels.Draw()
```
%% Output
%% Cell type:markdown id:defensive-doctrine tags:
## Generar toys a partir de un modelo
El modelo expuesto anteriormente sirve de base para generar "datos", que vamos a utilizar para el ajuste de los parámetros.
Los toys se generan utilizando un generador de números aleatorios de ROOT.
La siguiente celda genera entonces los toys de señal y fondo y los dibuja.
%% Cell type:code id:indonesian-charles tags:
``` python
nbins = 60
ROOT.gRandom.SetSeed(100) #importante fijar una semilla para la reproducibilidad de resultados.
toySignal = ROOT.gRandom.Poisson(nEntries*fSig) #el número de eventos señal y fondo lo obtenemos del valor de nuestro modelo y le permitimos que fluctue de forma Poissoniana.
toyBkg = ROOT.gRandom.Poisson(nEntries*(1-fSig))
print('Input signal events: {}'.format(toySignal))
histoBkg = ROOT.TH1D('histoBkg',"histoBkg",nbins, minVal,maxVal)
histoSig = ROOT.TH1D('histoSig',"histoSig",nbins, minVal,maxVal)
histoBkg.FillRandom('bkgModel',toyBkg)
histoSig.FillRandom('signalModel',toySignal)
cToys.cd(1)
histoSig.Draw('e')
cToys.cd(2)
histoBkg.Draw('e')
cToys.Draw()
```
%% Output
Input signal events: 537
%% Cell type:markdown id:stone-opening tags:
Y ahora, los combinamos y hacemos el ajuste de los parámetros.
%% Cell type:code id:coordinated-adjustment tags:
``` python
c1.cd()
#stack = ROOT.THStack('stack','stack')
#histoBkg.SetFillColor(2)
#histoSig.SetFillColor(4)
#stack.Add(histoBkg)
#stack.Add(histoSig)
pseudoData = ROOT.TH1D('pseudoData',"pseudoData",nbins, minVal,maxVal)
pseudoData.GetXaxis().SetTitle('x [GeV]')
pseudoData.GetYaxis().SetTitle('Entries / {} GeV'.format((maxVal-minVal)/nbins))
pseudoData.Add(histoBkg)
pseudoData.Add(histoSig)
pseudoData.Draw('e')
fullModel.SetParameters(nEntries*fSig,nEntries*(1-fSig),100.000000,125.000000,2.400000 ,50)
for i in range(6):
fullModel.ReleaseParameter(i)
fullModel.FixParameter(2,100)
pseudoData.Fit('fullModel')
#stack.Draw('same')
c1.Draw()
```
%% Output
FCN=63.0316 FROM MIGRAD STATUS=CONVERGED 136 CALLS 137 TOTAL
EDM=1.06449e-08 STRATEGY= 1 ERROR MATRIX ACCURATE
EXT PARAMETER STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 nsig 5.81424e+02 5.92347e+01 1.79365e-01 1.17869e-06
2 nbkg 9.62311e+03 1.83066e+02 3.98543e-01 -3.18049e-07
3 delta 1.00000e+02 fixed
4 mass 1.25194e+02 2.69660e-01 1.04322e-03 -2.35283e-04
5 sigma 2.54389e+00 2.77719e-01 8.94879e-04 -4.00985e-04
6 tau 4.68561e+01 1.39443e+00 3.13290e-03 -1.00685e-04
%% Cell type:markdown id:corporate-folks tags:
## ¿Qué estamos haciendo realmente?
En la parte anterior hemos obtenido una esimación de los parámetros y sus errores a través de una minimización del $\chi^{2}$ (en realidad no hemos sido nosotros, ROOT lo ha hecho utilizando un software que se llama MINUIT, y cuyo contenido sale un poco de la línea del curso..).
<h1><center>$\chi^{2} = \sum \frac{(x_{data}-x_{model})^{2}}{\sigma_{data}^{2}}$</center></h1>
Grosso modo: en el hiperespacio de parámetros libres (en nuestro caso, 5) ha generado una hipersuperficie de $\chi^{2}$. Lógicamente, cada conjunto de parámetros tiene un $\chi^{2}$
calculado el gradiente de la superficie en cada punto, y apuntado al punto en el cual es mínimo.
Para "verlo" mejor, vamos a hacer una simplificación: vamos a fijarnos exclusivamente en dos parámetros (la masa y la anchura de la señal) y vamos a ver qué cara tiene dicha superficie en ese subespacio de dos parámetros.
%% Cell type:raw id:aggregate-salad tags:
Primero creamos un histograma bidimensional
%% Cell type:code id:proof-mapping tags:
``` python
stepMass = 0.5
stepWidth = 0.25
binMass = int((maxVal-minVal)/stepMass)
binWidth = int(10./stepWidth)
chi2dPlot = ROOT.TH2D('chi2dPlot','chi2dPlot',binMass,minVal-stepMass/2.0,maxVal-stepMass/2.0,binWidth,0.5,10.5)
```
%% Cell type:markdown id:formed-particular tags:
Luego realizamos ajustes fijando tanto la masa como la anchura, de forma que nos permitan obtener el $\chi^{2}$ en función de dichos parámetros. Sin fijarlos, obtendríamos siempre el valor que más se ajuste a los datos (en el mínimo!).
%% Cell type:code id:liquid-rwanda tags:
``` python
#this is for the chi2 on the mass and width
chi2dPlot.Clear()
for masses in np.linspace(minVal, maxVal, num=binMass+1):
#print('Filling: mass and bin {} , {} '.format(masses,chi2dPlot.GetXaxis().FindBin(masses)))
for widths in np.linspace(0.25, 10.5, num=binWidth+1):
#print('Filling: mass and bin {} , {} '.format(widths,chi2dPlot.GetYaxis().FindBin(widths)))
fullModel.SetParameters(nEntries*fSig,nEntries*(1-fSig),100.000000,125.000000,2.400000 ,50)
fullModel.FixParameter(3,masses)
fullModel.FixParameter(4,widths)
pseudoData.Fit('fullModel','Q')
chi2dPlot.SetBinContent(chi2dPlot.FindBin(masses,widths),fullModel.GetChisquare())
```
%% Cell type:code id:patient-organ tags:
``` python
#remove borders..
chi2dPlot.SetAxisRange(0.5,10,"y")
chi2dPlot.SetAxisRange(111,150,"x")
chi2dPlot.GetXaxis().SetTitle('mass of the signal [GeV]')
chi2dPlot.GetYaxis().SetTitle('width of the signal [GeV]')
chi2dPlot.GetZaxis().SetTitle('#chi^{2}')
chi2dPlot.Draw('colz')
ROOT.gStyle.SetOptStat(000000)
c1.Draw()
```
%% Output
%% Cell type:code id:minute-blackberry tags:
``` python
```
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment