Skip to content
Snippets Groups Projects
Commit beb418ed authored by Francesco Santanastasio's avatar Francesco Santanastasio
Browse files

fit examples and counting experiment

parent fe8e30c2
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
CountingExperiment.ipynb 0 → 100644
+ 308
− 0
View file @ beb418ed
%% Cell type:code id:c6c01141 tags:
``` python
#plots will be shown inline
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
from scipy import stats
from scipy import special
from scipy import integrate
import pandas as pd
import math
import emcee
import corner
```
%% Cell type:markdown id:952245bd tags:
# Esperimento di conteggio in presenza di segnale e fondo con categorie
%% Cell type:code id:53316186 tags:
``` python
# test
nobs = np.array([0])
nbkg = np.array([0])
eff = np.array([1])
#0.7TeV
#inclusive
#nobs = np.array([14784,617])
#nbkg = np.array([14784,617])
#eff = np.array([0.266,0.124])
#tag+notag
#nobscat = np.array([3989,10795,39,578])
#nbkgcat = np.array([3989,10795,39,578])
#effcat = np.array([0.081,0.184,0.034,0.09])
#1TeV
#inclusive
#nobs = np.array([3984,90])
#nbkg = np.array([3984,90])
#eff = np.array([0.285,0.153])
#tag+notag
#nobscat = np.array([912,3072,6,84])
#nbkgcat = np.array([912,3072,6,84])
#effcat = np.array([0.093,0.192,0.051,0.102])
#2TeV
#inclusive
#nobs = np.array([67,0])
#nbkg = np.array([67,0])
#eff = np.array([0.316,0.2])
#tag+notag
#nobscat = np.array([25,42,0,0])
#nbkgcat = np.array([28,42,0,0])
#effcat = np.array([0.055,0.261,0.074,0.127])
#3TeV
#inclusive
#nobs = np.array([2,0])
#nbkg = np.array([2,0])
#eff = np.array([0.32,0.217])
#tag+notag
#nobscat = np.array([0,2,0,0])
#nbkgcat = np.array([0,2,0,0])
#effcat = np.array([0.02,0.3,0.033,0.184])
#4TeV
#inclusive
#nobs = np.array([0,0])
#nbkg = np.array([0.37,0])
#eff = np.array([0.309,0.2])
#tag+notag
#nobscat = np.array([0,0,0,0])
#nbkgcat = np.array([0.37,0,0,0])
#effcat = np.array([0.011,0.298,0.008,0.192])
#5TeV
#inclusive
#nobs = np.array([0,0])
#nbkg = np.array([0.1,0])
#eff = np.array([0.293,0.182])
#tag+notag
#nobscat = np.array([0,0,0,0])
#nbkgcat = np.array([0.1,0,0,0])
#effcat = np.array([0.015,0.278,0.012,0.170])
#lumi = 140 #fb-1
lumi =1 #fb-1
```
%% Cell type:code id:442b724f tags:
``` python
def log_likelihood(theta,nobs,nbkg,eff,lumi):
sigma = theta
nexp = eff * sigma * lumi + nbkg
#print (nexp)
return math.log(np.prod( stats.poisson.pmf(nobs,nexp) ))
```
%% Cell type:code id:1f8fd349 tags:
``` python
def log_prior(theta):
sigma = theta
if 0 < sigma < 10:
return 0.0
return -np.inf
```
%% Cell type:code id:1b76422e tags:
``` python
def log_probability(theta, nobs,nbkg,eff,lumi):
lp = log_prior(theta)
if not np.isfinite(lp):
return -np.inf
return lp + log_likelihood(theta, nobs,nbkg,eff,lumi)
```
%% Cell type:code id:a1c463f3 tags:
``` python
nwalkers = 500
niter = 500
initial = np.array([0.2])
ndim = len(initial)
p0 = [np.array(initial) + 1e-7 * np.random.randn(ndim) for i in range(nwalkers)]
#print (p0)
#no PPS categories
sampler = emcee.EnsembleSampler(nwalkers, ndim, log_probability, args=(nobs,nbkg,eff,lumi))
#with PPS categories
#sampler = emcee.EnsembleSampler(nwalkers, ndim, log_probability, args=(nobscat,nbkgcat,effcat,lumi))
sampler.run_mcmc(p0, niter, progress=True);
```
%% Output
100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 500/500 [00:15<00:00, 32.65it/s]
%% Cell type:code id:6bb9b203 tags:
``` python
flat_samples = sampler.get_chain(discard=100, thin=15, flat=True)
print(flat_samples.shape)
```
%% Output
(13000, 1)
%% Cell type:code id:b3bed1a2 tags:
``` python
labels = ["sigma"]
fig = corner.corner(
flat_samples, labels=labels
);
```
%% Output
%% Cell type:code id:d6babb25 tags:
``` python
from IPython.display import display, Math
for i in range(ndim):
mcmc = np.percentile(flat_samples[:, i], [16, 50, 84])
q = np.diff(mcmc)
txt = "\mathrm{{{3}}} = {0:.3f}_{{-{1:.3f}}}^{{{2:.3f}}}"
txt = txt.format(mcmc[1], q[0], q[1], labels[i])
display(Math(txt))
mcmc_up = np.percentile(flat_samples[:, i], [95])
txt_up = "\mathrm{{95\% \, CL \, Upper \, limit \, on \, {1}}} = {0:.3f} \mathrm{{\, fb}}"
txt_up = txt_up.format(mcmc_up[0], labels[i])
display(Math(txt_up))
```
%% Output
$\displaystyle \mathrm{sigma} = 0.666_{-0.497}^{1.155}$
$\displaystyle \mathrm{95\% \, CL \, Upper \, limit \, on \, sigma} = 3.017 \mathrm{\, fb}$
%% Cell type:code id:cfd9a67f tags:
``` python
```
LinerFitExamples.ipynb 0 → 100644
+ 680
− 0
View file @ beb418ed
This diff is collapsed.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment