Skip to content

Fixes to the log likelihood calculation #108

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
senkin merged 1 commit into from
Aug 1, 2014
Merged

Fixes to the log likelihood calculation #108

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
27 changes: 14 additions & 13 deletions tools/Fitting.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -219,8 +219,8 @@ def fit( self ):
# error flag for functions passed as reference.set to as 0 is no error
errorFlag = Long( 2 )

N_total = self.fit_data_collection.max_n_data() * 2
N_min = 0
N_max = self.fit_data_collection.max_n_data() * 2

param_index = 0

Expand All @@ -237,9 +237,9 @@ def fit( self ):
# >0 if MNPARM unable to implement definition
for sample in self.samples: # all samples but data
if self.n_distributions > 1:
gMinuit.mnparm( param_index, sample, self.normalisation[self.distributions[0]][sample], 10.0, N_min, N_total, errorFlag )
gMinuit.mnparm( param_index, sample, self.normalisation[self.distributions[0]][sample], 10.0, N_min, N_max, errorFlag )
else:
gMinuit.mnparm( param_index, sample, self.normalisation[sample], 10.0, N_min, N_total, errorFlag )
gMinuit.mnparm( param_index, sample, self.normalisation[sample], 10.0, N_min, N_max, errorFlag )
param_index += 1

arglist = array( 'd', 10 * [0.] )
Expand Down Expand Up @@ -306,18 +306,19 @@ def logLikelihood( self, nParameters, gin, f, par, iflag ):

def build_single_LL( self, data_vector, mc_vectors, normalisation, par ):
lnL = 0.0
for vector_entry, data_i in enumerate( data_vector ):
for vector_entry, v_data in enumerate( data_vector ):
x_i = 0
param_index = 0
for sample in self.samples:
x_i += par[param_index] * mc_vectors[sample][vector_entry]
self.param_indices[sample] = param_index
param_index += 1
data_term = normalisation[FitData.data_label] * data_i
if not data_term == 0 and not x_i == 0:
L = TMath.Poisson( data_term, x_i )
if not L == 0:
lnL += math.log( L )
data_i = normalisation[FitData.data_label] * v_data
L = TMath.Poisson( data_i, x_i )
if not L == 0:
lnL += math.log( L )
else:
lnL += -1e10
return lnL

def get_fit_normalisation_constraints( self, params ):
Expand Down Expand Up @@ -369,8 +370,8 @@ def fit( self ):
roofit_pdfs = {}
roofit_variables = {}

N_total = self.normalisation[self.data_label] * 2.
N_min = 0.
N_max = self.normalisation[self.data_label] * 2.
pdf_arglist = RooArgList()
variable_arglist = RooArgList()

Expand All @@ -386,7 +387,7 @@ def fit( self ):
roofit_variable = RooRealVar( sample, sample + " events",
self.normalisation[sample],
N_min,
N_total )
N_max )
roofit_variables[sample] = roofit_variable
pdf_arglist.add( roofit_pdf )
variable_arglist.add( roofit_variable )
Expand Down Expand Up @@ -476,13 +477,13 @@ def __init__( self, fit_data ):
samples = fit_data_1.samples
self.observables = {}
N_min = 0
N_total = fit_data_1.n_data() * 2
N_max = fit_data_1.n_data() * 2
for sample in samples:
self.observables[sample] = Observable( 'n_' + sample,
'number of ' + sample + " events",
fit_data_1.normalisation[sample],
N_min,
N_total,
N_max,
"events" )

# next create the models
Expand Down