def getPerformanceOnMC(
    files,
    histogramForEstimation="QCDStudy/PFIsolation_WithMETCutAndAsymJetCuts_DR03",
    bjetBin="",
    function="expo",
    fitRange=(0.3, 1.6),
    additionFitRanges=[(0.2, 1.6), (0.4, 1.6)],
):
    if bjetBin:
        histogramForEstimation = histogramForEstimation + "_" + bjetBin

    hists = [histogramForEstimation]
    hists = getHistsFromFiles(hists, files)
    hists = addSampleSum(hists)

    histogramForComparison = hists["qcd"][histogramForEstimation]
    histogramForEstimation = hists["allMC"][histogramForEstimation]

    estimate, absoluteError = getQCDEstimateFor(
        histogramForEstimation, function, fitRange=fitRange, additionFitRanges=additionFitRanges
    )

    qcdInSignalRegion, qcdError = getIntegral(histogramForComparison, (0, 0.1))

    N_est = ufloat((estimate, absoluteError))
    N_qcd = ufloat((qcdInSignalRegion, qcdError))

    relativeDeviation = N_est / N_qcd

    result = {}
    result["performance"] = (relativeDeviation.nominal_value, relativeDeviation.std_dev())
    result["estimate"] = (estimate, absoluteError)
    result["qcdInSignalRegion"] = (qcdInSignalRegion, qcdError)

    return result

def get_vhtt_yield(mass):
    log.warning("Get VHTT yield %s", mass)
    mmt_yield = results_file.Get('mmt_mumu_final_count/VH%i' % mass).Integral()
    emt_yield = results_file.Get('emt_emu_final_count/VH%i' % mass).Integral()

    _, (mmt_passed, mmt_total) = get_stat_error('VH%i' % mass, mmt_yield)
    _, (emt_passed, emt_total) = get_stat_error('VH%i' % mass, emt_yield)

    assert(mmt_total == emt_total)

    mmt_passed = ufloat( (mmt_passed, math.sqrt(mmt_passed) ))
    emt_passed = ufloat( (emt_passed, math.sqrt(emt_passed) ))
    total = ufloat( (mmt_total, math.sqrt(mmt_total) ))

    all_passed = mmt_passed + emt_passed

    log.warning("total %s", total)
    log.warning("all %s", all_passed)

    #multply by fraction of ZH_WH_ttH that is WH
    wh_xsec = ROOT.HiggsCSandWidth().HiggsCS(3, mass, 7)
    zh_xsec = ROOT.HiggsCSandWidth().HiggsCS(4, mass, 7)
    tth_xsec = ROOT.HiggsCSandWidth().HiggsCS(5, mass, 7)

    wh_fraction = wh_xsec/(wh_xsec + zh_xsec + tth_xsec)

    log.warning("wh_frac %s", wh_fraction)

    br_bonus = 1.0/(0.1075+0.1057+0.1125)

    log.warning("br_bonus %s", br_bonus)

    return 100*br_bonus*(all_passed/total)/wh_fraction

def split_dat(data):
    """
        Takes the row in the form r_n e_r_n ext_n e_ext_n and returns two lists
        one with the radii and the other with the extinction values both lists
        are of ufloats

        Parameters
        --------
                data : list
                A list of floats of the radii and errors and extinctions and 
                errors.

        Returns 
        --------
                rad : list
                List of ufloats of the radii.

                ext : list
                List of ufloats of the extinctions.
    """
    rad, ext = [], []

    for i in range(0, len(data) - 1, 4):
        vec = [data[i], data[i+1], data[i+2], data[i+3]]
        
        if vec[0].strip() != '':
            radius = ufloat(vec[0], vec[1])
            extinct = ufloat(vec[2], vec[3])

            rad.append(radius)
            ext.append(extinct)

    return ext, rad

 def fit_single_line_BS(self, x, y, zero_lev, err_continuum, fitting_parameters, bootstrap_iterations = 1000):
     
     #Declare parameters and containers for the fit
     params_dict     = fitting_parameters.valuesdict()
     initial_values  = [params_dict['A0'], params_dict['mu0'], params_dict['sigma0']]
     area_array      = empty(bootstrap_iterations)
     params_array    = empty([3, bootstrap_iterations])
     n_points        = len(x)
                         
     #Perform the fit
     for i in range(bootstrap_iterations):
         y_new               = y + np_normal_dist(0, err_continuum, n_points)
         area_array[i]       = simps(y_new, x) - simps(zero_lev, x)
         best_vals, covar    = curve_fit(gaussian_curveBS, (x, zero_lev), y_new, p0=initial_values, maxfev = 1600)
         params_array[:,i]   = best_vals
                 
     #Compute Bootstrap output
     mean_area, std_area = mean(area_array), std(area_array)
     mean_params_array, stdev_params_array = params_array.mean(1), params_array.std(1)
     
     #Store the data
     self.fit_dict['area_intg'],     self.fit_dict['area_intg_er']   = mean_area, std_area
     self.fit_dict['A0_norm'],       self.fit_dict['A0_norm_er']     = mean_params_array[0], stdev_params_array[0]
     self.fit_dict['mu0_norm'],      self.fit_dict['mu0_norm_er']    = mean_params_array[1], stdev_params_array[1]  
     self.fit_dict['sigma0_norm'],   self.fit_dict['sigma0_norm_er'] = mean_params_array[2], stdev_params_array[2]  
             
     A = ufloat(mean_params_array[0], stdev_params_array[0]) 
     sigma = ufloat(mean_params_array[2], stdev_params_array[2])         
     fwhm0_norm = 2.354820045 * sigma
     areaG0_norm = A * sigma * self.sqrt2pi
   
     self.fit_dict['fwhm0_norm'], self.fit_dict['fwhm0_norm_er'] = fwhm0_norm.nominal_value, fwhm0_norm.std_dev
     self.fit_dict['area_G0_norm'], self.fit_dict['area_G0_norm_er'] = areaG0_norm.nominal_value, areaG0_norm.std_dev
                        
     return

def line_fit(x, y, errors=True):
    """
    Simple helper function that speeds up line fitting. Uses lmfit
    
    Parameters
    ----------
    x, y (float)
        Sample length x, y arrays of data to be fitted

    Returns
    -------
    Returns slope and intercept, with uncertainties (use uncertainties package if availabe)
    Also returns fit object out, can be dropped
    """
    from lmfit.models import LinearModel
    mod = LinearModel()
    par = mod.guess(y, x=x)
    out = mod.fit(y, par, x=x)

    s = out.params['slope']
    i = out.params['intercept']

    if errors:
        try:
            from uncertainties import ufloat
            return ufloat(s.value, s.stderr), ufloat(i.value, i.stderr), out
        except:
            return s.value, s.stderr, i.value, i.stderr, out
    else:
        return s.value, i.value, out

def calculate_latencies(version_dates):
    linux_latencies = latency(version_dates['linux'], OrderedDict(avo.os_to_kernel))
    set_latex_value('linuxMeanUpdateLatency', ufloat(statistics.mean(linux_latencies.values()),statistics.stdev(linux_latencies.values())))
    openssl_latencies = latency(version_dates['openssl'], OrderedDict(avo.os_to_project['openssl']))
    set_latex_value('opensslMeanUpdateLatency', ufloat(statistics.mean(openssl_latencies.values()),statistics.stdev(openssl_latencies.values())))
    bouncycastle_latencies = latency(version_dates['bouncycastle'], OrderedDict(avo.os_to_project['bouncycastle']))
    set_latex_value('bouncycastleMeanUpdateLatency',ufloat(statistics.mean(bouncycastle_latencies.values()),statistics.stdev(bouncycastle_latencies.values())))

def combine_complex_df( df1, df2 ):
	'''
	Takes a 2 pandii dataframes of the form:
		A 	| 	B   	 	A 	| 	B 
	  (v,e) | (v,e)		  (v,e) | (v,e) 

	Returns 1 pandas dataframe of the form
			      A   |   B 
	  			(v,e) | (v,e)
	'''
	from uncertainties import ufloat
	l1=df1.columns.tolist()
	l2=df2.columns.tolist()
	if l1 != l2:
		print "Trying to combine two non compatible dataframes"
		print l1
		print l2
		return

	combined_result = {}
	for sample in l1:
		results = []
		for entry1, entry2 in zip(df1[sample], df2[sample]):
			v1 = ufloat(entry1[0], entry1[1])
			v2 = ufloat(entry2[0], entry2[1])
			s = v1 + v2
			results.append( ( s.nominal_value, s.std_dev ) )
		combined_result[sample] = results
	df = dict_to_df(combined_result)
	return df

def calculate_flux(f, age, arar_constants=None):
    """
        #rad40: radiogenic 40Ar
        #k39: 39Ar from potassium
        f: F value rad40Ar/39Ar
        age: age of monitor in years

        solve age equation for J
    """
    # if isinstance(rad40, (list, tuple)):
    # rad40 = ufloat(*rad40)
    # if isinstance(k39, (list, tuple)):
    # k39 = ufloat(*k39)

    if isinstance(f, (list, tuple)):
        f = ufloat(*f)

    if isinstance(age, (list, tuple)):
        age = ufloat(*age)
        #    age = (1 / constants.lambdak) * umath.log(1 + JR)
    try:
        # r = rad40 / k39
        if arar_constants is None:
            arar_constants = ArArConstants()

        j = (umath.exp(age * arar_constants.lambda_k.nominal_value) - 1) / f
        return j.nominal_value, j.std_dev
    except ZeroDivisionError:
        return 1, 0

def small_sep02 ( mx , mxErr, lag=0):
	""" 
	Given a frequency matrix and errors calculates the (scaled) small 
	separation d02 and propagates errors.
	
	Notes
	-----
	The parameter lag is the difference between the radial orders of
	the first modes of l=0 and l=2.  
	"""

	(Nn, Nl) = np.shape(mx)
	d02    = np.zeros((1,Nn))
	d02Err = np.zeros((1,Nn))
	d02.fill(None)      # values that can't be calculated are NaN

	
	for n in range(1-lag,Nn):
		if (mx[n,0] != 0. and mx[n-1+lag,2] != 0.):
		    a = un.ufloat( (mx[n,0], mxErr[n,0]) )
		    b = un.ufloat( (mx[n-1+lag,2], mxErr[n-1+lag,2]) )
		    result = (a-b) / 3.	
		    d02[0,n-1+lag]    = un.nominal_value(result)
		    d02Err[0,n-1+lag] = un.std_dev(result)
		
	return d02, d02Err

def day3_vacuum():
    plt.clf()
    ux = unp.uarray([800, 1000, 500], ux_error)  # V
    ug_crit = unp.uarray([110, 140, 30], ug_error)  # V
    omega = 2 * math.pi * ufloat(48, 1)  # Hz

    ug_crit = _apply_additional_proportional_error(ug_crit, stability_uncertainty_vacuum)
    ux *= ux_correction
    ug_crit *= ug_correction

    stability(ux, ug_crit, omega, label="Particle V1")

    plt.title("p = 300 mbar")
    plt.legend(loc=2)
    plt.savefig("images/stability_vacuum_1.pdf")

    plt.clf()
    ux = unp.uarray([700, 800], ux_error)  # V
    ug_crit = unp.uarray([100, 140], ug_error)  # V
    omega = 2 * math.pi * ufloat(45, 1)  # Hz

    ug_crit = _apply_additional_proportional_error(ug_crit, stability_uncertainty_vacuum)
    ux *= ux_correction
    ug_crit *= ug_correction

    stability(ux, ug_crit, omega, label="Particle V2")

    plt.title("p = 180 mbar")
    plt.legend(loc=2)
    plt.savefig("images/stability_vacuum_2.pdf")

def safe_div(n, d):
    if n.size == 0:
        return None
    if n.shape != d.shape:
        return None
    if isinstance(n[0, 0], UFloat) is False:
        return None
    z = np.zeros_like(n)
    z.fill(ufloat(np.nan, np.nan))
    it = np.nditer(d, op_flags=["readonly"], flags=["multi_index", "refs_ok"])
    while not it.finished:
        i = it.multi_index[0]
        j = it.multi_index[1]
        if abs(d[i, j].n) == 0.0:
            # z[i,j] = ufloat(0.,0.)
            z[i, j] = ufloat(np.nan, np.nan)
        #            if abs(n[i,j].n) > 0. :
        #                z[i,j] = ufloat(np.nan,np.nan)
        #            else :
        #                err = np.sqrt( n[i,j].s**2. + d[i,j].s**2. )
        #                z[i,j] = ufloat(0.,err)
        # z[i,j] = ufloat(np.nan,np.nan)
        # z[i,j] = ufloat(np.inf,np.inf)
        else:
            z[i, j] = n[i, j] / d[i, j]
        it.iternext()
    return z

def day2_air():
    plt.clf()
    ux = unp.uarray([1000, 860, 800, 760, 1000, 920, 620], ux_error)  # V
    ug_crit_1 = unp.uarray([390, 170, 600, 120, 225, 175, 53], ug_error)  # V
    ug_crit_2 = unp.uarray([390, 170, 600, 120, 160, 145, 67], ug_error)  # V
    ug_crit_3 = unp.uarray([390, 370, 600, 120, 225, 175, 67], ug_error)  # V
    omega = 2 * math.pi * ufloat(28, 1)  # Hz

    ux *= ux_correction
    for i, ug_crit in enumerate((ug_crit_1, ug_crit_2, ug_crit_3), 1):
        ug_crit = _apply_additional_proportional_error(ug_crit, stability_uncertainty_air)
        ug_crit *= ug_correction

        stability(ux, ug_crit, omega, "Particle A%d" % i)

    plt.legend(loc=2)
    plt.savefig("images/stability_air_1.pdf")

    plt.clf()
    ux = unp.uarray([550, 740, 870, 1040], ux_error)
    ug_crit = unp.uarray([39, 84, 133, 147], ug_error)
    omega = 2 * math.pi * ufloat(32, 1)  # Hz

    ug_crit = _apply_additional_proportional_error(ug_crit, stability_uncertainty_air)
    ux *= ux_correction
    ug_crit *= ug_correction

    stability(ux, ug_crit, omega, label="Particle A4")

    plt.legend(loc=2)
    plt.savefig("images/stability_air_2.pdf")

def get_1d_loose_efficiencies( int_stat, int_syst, lead_reg, lead_ptrange, systematics=None) :

    eff_stat = {}
    eff_syst = {}

    if int_stat['lead']['real']['loose'].n == 0 :
        eff_stat['eff_R_T_lead'] = ufloat( 1.0, int_stat['lead']['real']['tight'].s/int_stat['lead']['real']['tight'].n )
    else :
        eff_stat['eff_R_T_lead'] = int_stat['lead']['real']['tight'] / (int_stat['lead']['real']['tight']+int_stat['lead']['real']['loose'])
    eff_stat['eff_F_T_lead'] = int_stat['lead']['fake']['tight'] / (int_stat['lead']['fake']['tight']+int_stat['lead']['fake']['loose'])
    eff_stat['eff_R_L_lead'] = int_stat['lead']['real']['loose'] / (int_stat['lead']['real']['tight']+int_stat['lead']['real']['loose'])
    eff_stat['eff_F_L_lead'] = int_stat['lead']['fake']['loose'] / (int_stat['lead']['fake']['tight']+int_stat['lead']['fake']['loose'])

    eff_syst['eff_R_T_lead'] = int_syst['lead']['real']['tight'] / (int_syst['lead']['real']['tight']+int_syst['lead']['real']['loose'])
    eff_syst['eff_F_T_lead'] = int_syst['lead']['fake']['tight'] / (int_syst['lead']['fake']['tight']+int_syst['lead']['fake']['loose'])
    eff_syst['eff_R_L_lead'] = int_syst['lead']['real']['loose'] / (int_syst['lead']['real']['tight']+int_syst['lead']['real']['loose'])
    eff_syst['eff_F_L_lead'] = int_syst['lead']['fake']['loose'] / (int_syst['lead']['fake']['tight']+int_syst['lead']['fake']['loose'])


    # Do systematics
    # the integrals may already have some systematics
    # that are propagated to the eff_*
    # therefore, don't overwrite the existing 
    # systematics, but make a ufloat with a
    # zero value, and non-zero syst
    eff_syst['eff_R_L_lead'] = eff_syst['eff_R_L_lead'] + ufloat( 0.0, math.fabs(eff_syst['eff_R_L_lead'].n)*get_syst_uncertainty( 'RealTemplateNom', lead_reg, lead_ptrange, 'real', 'loose' ), 'Template_lead_real_loose')
    eff_syst['eff_F_L_lead'] = eff_syst['eff_F_L_lead'] + ufloat( 0.0, math.fabs(eff_syst['eff_F_L_lead'].n)*get_syst_uncertainty( 'FakeTemplate%s'%systematics, lead_reg, lead_ptrange, 'fake', 'loose' ), 'Template_lead_fake_loose' )
    eff_syst['eff_R_T_lead'] = eff_syst['eff_R_T_lead'] + ufloat( 0.0, math.fabs(eff_syst['eff_R_T_lead'].n)*get_syst_uncertainty( 'RealTemplateNom', lead_reg, lead_ptrange, 'real', 'loose' ), 'Template_lead_real_loose')
    eff_syst['eff_F_T_lead'] = eff_syst['eff_F_T_lead'] + ufloat( 0.0, math.fabs(eff_syst['eff_F_T_lead'].n)*get_syst_uncertainty( 'FakeTemplate%s'%systematics, lead_reg, lead_ptrange, 'fake', 'loose' ), 'Template_lead_fake_loose' )

    return eff_stat, eff_syst

    def _set_age_values(self, f, include_decay_error=False):
        arc = self.arar_constants
        j = copy(self.j)
        if j is None:
            j = ufloat(1e-4, 1e-7)
        j.tag = 'Position'
        j.std_dev = self.position_jerr or 0
        age = age_equation(j, f, include_decay_error=include_decay_error, arar_constants=arc)
        self.uage_w_position_err = age

        j = self.j
        if j is None:
            j = ufloat(1e-4, 1e-7, tag='J')

        age = age_equation(j, f, include_decay_error=include_decay_error, arar_constants=arc)
        self.uage_w_j_err = age

        j = copy(self.j)
        if j is None:
            j = ufloat(1e-4, 1e-7, tag='J')

        j.std_dev = 0
        age = age_equation(j, f, include_decay_error=include_decay_error, arar_constants=arc)
        self.uage = age

        self.age = nominal_value(age)
        self.age_err = std_dev(age)
        self.age_err_wo_j = std_dev(age)

        for iso in self.itervalues():
            iso.age_error_component = self.get_error_component(iso.name)

 def test_decimal_places(self):
     cases = [(0, '0'), (-1, '-1'), (0.012, '0.01'), (0.016, '0.02')
             , (ufloat(0,0), r'$0.0 \pm 0.0$'), (ufloat(10,0), r'$10.0 \pm 0.0$')
             , (ufloat(0.12, 0.012), r'$0.12 \pm 0.01$'), (ufloat(0.12, 0.016), r'$0.12 \pm 0.02$')
             , (ufloat(-0.12, 0.016), r'$-0.12 \pm 0.02$')]
     for test_case, result in cases:
         self.assertEqual(display_num(test_case, decimal_places=2), result)

def small_sep13 ( mx , mxErr, lag=0):
	""" 
	Given a frequency matrix and errors calculates the (scaled) small 
	separation d13 and propagates errors.
	
	Notes
	-----
	The parameter lag is the difference between the radial orders of
	the first modes of l=1 and l=3.  
	"""

	(Nn, Nl) = np.shape(mx)
	d13    = np.zeros((1,Nn))
	d13Err = np.zeros((1,Nn))
	d13.fill(None)      # values that can't be calculated remain NaN

	for n in range(1-lag,Nn):
		if (mx[n,1] != 0. and mx[n-1+lag,3] != 0.):     
        		a = un.ufloat( (mx[n,1], mxErr[n,1]) )
        		b = un.ufloat( (mx[n-1+lag,3], mxErr[n-1+lag,3]) )
        		result = (a-b) / 5.
        		d13[0,n-1+lag]    = un.nominal_value(result)
        		d13Err[0,n-1+lag] = un.std_dev(result)

	return d13, d13Err

    def eval(self, vertex_generator, nevals, nreps=16, ndaq=50):
        """
        Return the negative log likelihood that the detector event set in the
        constructor or by set_event() was the result of a particle generated
        by `vertex_generator`.  If `nreps` set to > 1, each set of photon
        vertices will be propagated `nreps` times.
        """
        ntotal = nevals * nreps * ndaq

        hit_prob, pdf_prob, pdf_prob_uncert = \
            self.eval_channel_vbin(vertex_generator, nevals, nreps, ndaq)

        # NLL calculation: note that negation is at the end
        # Start with the probabilties of hitting (or not) the channels
        
        # Flip probability for channels in event that were not hit
        hit_prob[~self.event.channels.hit] = 1.0 - hit_prob[~self.event.channels.hit]
        # Apply a floor so that we don't take the log of zero
        hit_prob = np.maximum(hit_prob, 0.5 / ntotal)

        hit_channel_prob = np.log(hit_prob).sum()
        log_likelihood = ufloat((hit_channel_prob, 0.0))

        # Then include the probability densities of the observed
        # charges and times.
        log_likelihood += ufloat((np.log(pdf_prob[self.event.channels.hit]).sum(),
                                  0.0))
        
        return -log_likelihood

    def get_mean_raw(self, tau=None):
        vs = []
        corrfunc = self._deadtime_correct
        for r in six.itervalues(self._cp):
            n = int(r['NShots'])
            nv = ufloat(float(r['Ar40']), float(r['Ar40err'])) * 6240
            dv = ufloat(float(r['Ar36']), float(r['Ar36err'])) * 6240
            if tau:
                dv = corrfunc(dv, tau * 1e-9)

            vs.append((n, nv / dv))

        key = lambda x: x[0]
        vs = sorted(vs, key=key)

        mxs = []
        mys = []
        mes = []
        for n, gi in groupby(vs, key=key):
            mxs.append(n)
            ys, es = list(zip(*[(nominal_value(xi[1]), std_dev(xi[1])) for xi in gi]))

            wm, werr = calculate_weighted_mean(ys, es)
            mys.append(wm)
            mes.append(werr)

        return mxs, mys, mes

def interference_corrections(a39, a37,
                             production_ratios,
                             arar_constants=None,
                             fixed_k3739=False):
    if production_ratios is None:
        production_ratios = {}

    if arar_constants is None:
        arar_constants = ArArConstants()

    pr = production_ratios
    k37 = ufloat(0, 0, tag='k37')

    if arar_constants.k3739_mode.lower() == 'normal' and not fixed_k3739:
        # iteratively calculate 37, 39
        for _ in range(5):
            ca37 = a37 - k37
            ca39 = pr.get('Ca3937', 0) * ca37
            k39 = a39 - ca39
            k37 = pr.get('K3739', 0) * k39
    else:
        if not fixed_k3739:
            fixed_k3739 = arar_constants.fixed_k3739

        ca37, ca39, k37, k39 = apply_fixed_k3739(a39, pr, fixed_k3739)

    k38 = pr.get('K3839', 0) * k39

    if not arar_constants.allow_negative_ca_correction:
        ca37 = max(ufloat(0, 0), ca37)

    ca36 = pr.get('Ca3637', 0) * ca37
    ca38 = pr.get('Ca3837', 0) * ca37

    return k37, k38, k39, ca36, ca37, ca38, ca39

        def func(ra):
            u, l = ra.split(':')
            try:
                ru = self.signals[u]
                rl = self.signals[l]
            except KeyError:
                return ''

            if cfb:
                bu = ufloat(0, 0)
                bl = ufloat(0, 0)
                try:
                    bu = self.baselines[u]
                    bl = self.baselines[l]
                except KeyError:
                    pass
                try:
                    rr = (ru - bu) / (rl - bl)
                except ZeroDivisionError:
                    rr = ufloat(0, 0)
            else:
                try:
                    rr = ru / rl
                except ZeroDivisionError:
                    rr = ufloat(0, 0)

            res = '{}/{}={} '.format(u, l, pad('{:0.4f}'.format(rr.nominal_value))) + \
                  PLUSMINUS + pad(format('{:0.4f}'.format(rr.std_dev)), n=6) + \
                    self._get_pee(rr)
            return res