A short minimal example showcasing an issue with RooFit. You need ROOT with pyroot to run it.

roofit_problem

import scipy.integrate
import numpy
import matplotlib.pyplot as plt

points = [ 4425.941,  4369.934,  5903.859,  4432.979,  6035.395,  4685.236,  4382.774, 4532.092,  4591.621,  4629.322,  4739.585,  4529.45,   4500.441,  5752.579]
ranges = [(4367.,5020.),(5626,6366)]


def plot(exponent):
	_, borders, _ = plt.hist(points)
	binwidth = borders[1] - borders[0] #avoid potential trouble with border of a possible underflow bin
	x = numpy.linspace(ranges[0][0], ranges[1][1])
	plt.plot(x, numpy.exp(exponent*x)/get_norm(exponent) * binwidth * len(points))
	plt.show()

def get_norm(exponent):
	norm = 0
	for range in ranges:
		integral, error = scipy.integrate.quad(lambda x: numpy.exp(exponent * x), range[0], range[1])
		norm += integral
	return norm


def fit_root():
	import ROOT
	mass = ROOT.RooRealVar("B_s0_MM","B_s0_MM", ranges[0][0], ranges[1][1])
	mass.setRange("fitRangeLeft",  ranges[0][0] , ranges[0][1])
	mass.setRange("fitRangeRight", ranges[1][0], ranges[1][1])

	exponent = ROOT.RooRealVar("exponent","exponent", -0.0014, -0.1, 0.)
	model = ROOT.RooExponential("data_model", "data_model", mass, exponent)

	data = ROOT.RooDataSet("data","data", ROOT.RooArgSet(mass))
	for point in points:
		mass.setVal(point)
		data.addFast(ROOT.RooArgSet(mass))

	rooFitResult = model.fitTo(data, ROOT.RooFit.Save(), ROOT.RooFit.Range("fitRangeLeft,fitRangeRight"))
	rooFitResult.Print()

def get_llh_value(exponent):
	norm = get_norm(exponent)
	llh = 1
	for point in points:
		llh *= numpy.exp(exponent * point)/norm

	print "LLH value for exponent " + str(exponent)+ " : " + str(llh)


#######################################
#
# The issue (TM)
#
#######################################

# first, fit the data with root
fit_root() 
# this returns a best fit exponent of -0.0045887

# calculate the likelihood value at roots best fit value:
get_llh_value(exponent = -0.0045887) 
# it gives a likelihood of 7.16715257847e-46

# compare to the best fit value at -0.001497 (which is the actual optimum)
get_llh_value(exponent = -0.001497 ) 
# it gives a likelihood of 1.82358886152e-42
# the likelihood is higher!

# This means root is not finding the optimum

plot(exponent = -0.0045887)
# It may also be interesting to look at the plot of roots best fit