def PlotStep(self, Primal,Primal0,col,dt = 0.2):
     
     Nshooting = self.Nshooting
     Nturbine = self.Nturbine
     
     time = {'States': [dt*k for k in range(Nshooting+1)],
             'Inputs': [dt*k for k in range(Nshooting)]}
 
     Nsubp = np.ceil(np.sqrt(Nturbine))
     key_dic = {'States': ['Og','beta'], 'Inputs': ['Tg']}
     counter = 0
     for type_ in key_dic.keys():
         for key in key_dic[type_]:
             for k in range(Nturbine):
                 Primal_struc  = self._TurbineV(Primal['Turbine'][k])
                 Primal0_struc = self._TurbineV(Primal0['Turbine'][k])
                 diff   = veccat(Primal_struc[type_,:,key])-veccat(Primal0_struc[type_,:,key])
            
                 plt.figure(10+counter)
                 plt.subplot(Nsubp,Nsubp,k)
                 plt.hold('on')
                 if (type_ == 'States'):
                     plt.plot(time[type_],diff,color=col)
                 else:
                     plt.step(time[type_],diff,color=col)    
         
             plt.title('Step'+key)   
             counter += 1

def mask_spectrum(flux_to_fit,interactive=True,mask_lower_limit=None,mask_upper_limit=None):

    """
    Interactively and iteratively creates a Boolean mask for a spectrum.

    """
    if interactive:
      plt.ion()
      continue_parameter = 'no'
      mask_switch = 'yes'
      while mask_switch == 'yes':
        pixel_array = np.arange(len(flux_to_fit))
        plt.figure()
        plt.step(pixel_array,flux_to_fit)
        mask_lower_limit_string = raw_input("Enter mask lower limit (in pixels): ")
        mask_lower_limit = float(mask_lower_limit_string)
        mask_upper_limit_string = raw_input("Enter mask upper limit (in pixels): ")
        mask_upper_limit = float(mask_upper_limit_string)
        mask = (pixel_array >= mask_lower_limit) & (pixel_array <= mask_upper_limit)
        flux_to_fit_masked = np.ma.masked_where(mask,flux_to_fit)
        plt.step(pixel_array,flux_to_fit_masked)
        continue_parameter = raw_input("Happy? (yes or no]): ")
        plt.close()

        if continue_parameter == 'yes':
          mask_switch = 'no'

    else:
      pixel_array = np.arange(len(flux_to_fit))
      mask = (pixel_array >= mask_lower_limit) & (pixel_array <= mask_upper_limit)

    return mask

def plotSegments(jason, name="", labels=False):
    # read delimeters for each stop
    delims = readDelims()

    stopAdds = [getStopSet(jason["packets"], stop) for stop in delims]
    intersectBins = [set() for stop in delims]
    combobs = combinations(range(len(delims)), 2)
    for combob in combobs:
        curSect = stopAdds[combob[0]].intersection(stopAdds[combob[1]])
        for i in range(combob[0], combob[1]):
            intersectBins[i+1].update(curSect)
    probably_junk = stopAdds[0].intersection(stopAdds[-1])

    y = [len(bin - probably_junk) + 2 for bin in intersectBins]
    y[0] = y[1]  # For labelling purposes

    x = [stop['start'] for stop in delims]
    realy = [stop['actual'] for stop in delims]
    plot.xlabel('Seconds since '+jason["initial_time"])
    plot.ylabel('Number of bus occupants (predicted)')
    plot.xlim(0, delims[-1]['end'])
    plot.title(name)
    plot.step(x, y)
    plot.step(x, realy, color="purple", where="post")

    if(labels):
        for stop in delims:
            annotate(stop["code"], stop["start"], stop["actual"], 10, 10)

    makeWidePlot("bus", "segments")

    plot.show()

def draw_eta_profile(analyzed_root_files, # path to edm-analyzed root files
                    gen_energy, # energy with which particle is generated
                    mylabel = "", # label for histograms
                    histcolor = '#000000'):
    
######################################################################################
# Function reads in the root files with edm-analyzed data,                           #
# chains the data from the different files to one root tree,                         #
# reads the relevant data to arrays                                                  #
# and creates and draws histograms.                                                  #
# Drawing of errobars is done here.                                                  #
# Settings for decorations such as axeslabels are not done here, but in main program #
######################################################################################
    t0 = time.time()
    ## add to tree all edm-analyzed reco files
    rh_tree = ROOT.TChain("demo/rh_tree")
    rh_tree.Add(analyzed_root_files)
    print "=========================================================="
    print "Events in",str(analyzed_root_files),": ", rh_tree.GetEntries()
    
    ## initialize arrays
    particle_etas = [] # list with all eta values in cut intervall
    eta_energies = [] # list with all corresponding total RecHit energies for each eta over all events
    
    ## loop over all gen_particles
    for i in range(rh_tree.GetEntries()):
        rh_tree.GetEntry(i)
        momentum_vector = ROOT.TVector3(rh_tree.gen_part_momentum_x,
                                        rh_tree.gen_part_momentum_y,
                                        rh_tree.gen_part_momentum_z)
        eta = momentum_vector.Eta()
        ## do eta cut

        if ((eta > min_eta) and (eta < max_eta)):
            # and ((rh_tree.ecal_total_energy+rh_tree.hcal_total_energy) < 40.)):
            particle_etas.append(eta)
            eta_energies.append(rh_tree.ecal_total_energy
                                + rh_tree.hcal_total_energy)

    particle_etas = np.array(particle_etas,dtype=float)
    eta_energies = np.array(eta_energies,dtype=float)
    cut_event_numb = particle_etas.size
    

    print "time needed to read in data and loop over entries: ", np.round(time.time()-t0,2)
    ## draw histogram
    # etas are histogrammed automatically, corresponding energies are used as weights

    binvalues, binedges, binerrors = calc_histogram(cut_event_numb,
                                                    particle_etas,
                                                    eta_energies/float(gen_energy))
    t1 = time.time()
    plt.step(binedges, np.append(binvalues, binvalues[-1]),
             where="post", color=histcolor,
             label=str(mylabel))
    plt.errorbar(np.array(binedges[:-1])+(binedges[1]-binedges[0])/2, binvalues,
                yerr=binerrors,
                fmt='none', ecolor=histcolor)
    print "time needed to draw plot: ", np.round(time.time()-t1,2)
    return binvalues, binerrors, binedges

    def testExponManyEvents(self):
        """
        generate and fit an exponential distribution with lifetime of 25
        make a plot in testExponManyEvents.png
        """
        tau = 25.0
        nBins = 400
        size = 100
        taulist = []
        for i in range(1000):
            x = range(nBins)
            timeHgValues = np.zeros(nBins, dtype=np.int64)
            timeStamps = expon.rvs(loc=0, scale=tau, size=size)
            ts64 = timeStamps.astype(np.uint64)
            tsBinner.tsBinner(ts64, timeHgValues)
            
            param = expon.fit(timeStamps)
            fit = expon.pdf(x,loc=param[0],scale=param[1])
            fit *= size
            print "i=",i," param[1]=",param[1]
            taulist.append(param[1]) 

        hist,bins = np.histogram(taulist, bins=20, range=(15,25))
        width = 0.7*(bins[1]-bins[0])
        center = (bins[:-1]+bins[1:])/2
        plt.step(center, hist, where = 'post')
        plt.savefig(inspect.stack()[0][3]+".png")

def values_fromCDF():
    '''Calculate an empirical cumulative distribution function, compare it with the exact one, and
    find the exact point for a specific data value.'''
    
    # Generate normally distributed random data
    myMean = 5
    mySD = 2
    numData = 100
    data = stats.norm.rvs(myMean, mySD, size=numData)
    
    # Calculate the cumulative distribution function, CDF
    numbins = 20
    counts, bin_edges = np.histogram(data, bins=numbins, normed=True)
    cdf = np.cumsum(counts)
    cdf /= np.max(cdf)
    
    # compare with the exact CDF
    plt.step(bin_edges[1:],cdf)
    plt.hold(True)
    plt.plot(x, stats.norm.cdf(x, myMean, mySD),'r')
    
    # Find out the value corresponding to the x-th percentile: the
    # "cumulative distribution function"
    value = 2
    myMean = 5
    mySD = 2
    cdf = stats.norm.cdf(value, myMean, mySD)
    print(('With a threshold of {0:4.2f}, you get {1}% of the data'.format(value, round(cdf*100))))
    
    # For the percentile corresponding to a certain value: 
    # the "inverse cumulative distribution function" 
    value = 0.025
    icdf = stats.norm.isf(value, myMean, mySD)
    print(('To get {0}% of the data, you need a threshold of {1:4.2f}.'.format((1-value)*100, icdf)))
    plt.show()

def main():
    '''The data in this example give the life talbe for motion sickness data
    from an experiment with vertical movement at a frequency of 0.167 Hz and
    acceleration 0.111 g, and of a second experiment with 0.333 Hz and acceleration
    of 0.222 g.
    '''
    
    # get the data
    data1 = getData('altman_13_2.txt', subDir='..\Data\data_altman')
    data2 = getData('altman_13_3.txt', subDir='..\Data\data_altman')
    
    # Determine the Kaplan-Meier curves
    (p1, r1, t1, sp1,se1) = kaplanmeier(data1)
    (p2, r2, t2, sp2,se2) = kaplanmeier(data2)
    
    # Make a combined plot for both datasets
    plt.step(t1,sp1, where='post')
    plt.hold(True)
    plt.step(t2,sp2,'r', where='post')
    
    plt.legend(['Data1', 'Data2'])
    plt.ylim(0,1)
    plt.xlabel('Time')
    plt.ylabel('Survival Probability')
    plt.show()
    
    # Check the hypothesis that the two survival curves are the same
    # --- >>> START stats <<< ---
    (p, X2) = logrank(data1, data2)
    # --- >>> STOP stats <<< ---
    
    return p    # supposed to be 0.073326322306832212

def fvc_plot_setup(hist_data, hist, binEdges, xlabel, title = ""):
    '''
    Plot the histogram, with removed observations highlighted
    
    :param array hist_data: raw values which have been binned to create hist
    :param array hist: values of histogram
    :param array binEdges: location of LH bin edge
    :param str xlabel: label for x-axis
    :param str title: title of plot
    
    :returns:
        plot-hist - useful histogram data to plot in log-scale
        bincenters - locations of centres of bins
    '''
    import matplotlib.pyplot as plt

    plot_hist = np.array([float(x) if x != 0 else 1e-1 for x in hist])
    plt.clf()
    bincenters = 0.5 * (binEdges[1:] + binEdges[:-1])
    plt.step(bincenters, plot_hist, 'b-', label = 'observations', where='mid')
            
    fit = utils.fit_gaussian(bincenters, hist, max(hist), mu = np.mean(hist_data), sig = np.std(hist_data))
    plot_gaussian = utils.gaussian(bincenters, fit)
    plt.plot(bincenters, plot_gaussian, 'r-', label = 'Gaussian fit')
    # sort labels and prettify
    plt.xlabel(xlabel)                    
    plt.ylabel("Frequency")
    plt.gca().set_yscale('log')
    plt.ylim([0.1,10000])
    plt.title(title)

    return plot_hist, bincenters # fvc_plot_setup

def dopsd(nfft = None, rpt = 10, plotdB=True):
  """
  Takes a snapshot, then computes, plots and writes out the Power Spectral
  Density functions.  The psd function is written into a file named "psd".
  This file will be overwritten with each call.  Arguments:

  nfft The number of points in the psd function.  Defaults to the number of
       points in a snapshot, the maximum which should be used.

  rpt  The numper of mesurements to be averaged for the plot and output file. 
       Defaults to 10.

  plotdB controls whether the plot is linear in power or in dB
  """
  if nfft == None:
    nfft = numpoints
  for i in range(rpt):
    power, freqs = adc5g.get_psd(roach2, snap_name, samp_freq*1e6, 8, nfft)
    if i == 0:
      sp = power
    else:
      sp += power
  sp /= rpt
  if plotdB:
    plt.step(freqs, 10*np.log10(sp))
  else:
    plt.step(freqs, (sp))
  plt.show(block = False)
  data = np.column_stack((freqs/1e6, 10*np.log10(sp)))
  np.savetxt("psd", data, fmt=('%7.2f', '%6.1f'))

def plot_pdf(n):
	w = 0.5
	r = 0.9
	X = np.random.uniform(0,1,n)
	Y = map(lambda x: w*x + (1-w)*r, X)
	# CDF estimate of X
	plt.figure()
	ecdf_x = sm.tools.ECDF(X)
	x = np.linspace(min(X), max(X))
	f_x = ecdf_x(x)
	plt.step(x, f_x)
	plt.xlabel(r"$x$")
	plt.ylabel(r"Empirical CDF of $X \sim U(0,1)$")
	plt.grid()
	# PDF estimate of Y
	plt.figure()
	count, bins, patches = plt.hist(Y, n/200, normed=1)
	# CDF estimate of Y
	plt.figure()
	ecdf_y = sm.tools.ECDF(Y)
	y = np.linspace(min(Y), max(Y))
	f_y = ecdf_y(y)
	plt.step(y, f_y)
	plt.xlabel(r"$y$")
	plt.ylabel(r"Empirical CDF of $Y=w\cdot X + (1-w)\cdot r$")
	plt.grid()

def pltLumFun(data,lumbins,color='blue',linestyle='-',redshift=1,overplot=False,plotdata=True,label=None,linewidth=2):
	zz, = np.where(z==redshift)
	plt.step(lumbins,np.log10(np.append(data[zz,:],data[zz,-1])),color=color,linestyle=linestyle,label=label,lw=linewidth)
	if plotdata==True:
		obs = readcol.readcol('/nobackupp8/mtremmel/DATA/QSOdata/bol_lf_point_dump.dat',twod=False,asdict=True,skipline=38)
                obs2 = readcol.readcol('/nobackupp8/mtremmel/DATA/QSOdata/M1450z5_McGreer13.dat',twod=False,asdict=True,skipline=1)
		tt, = np.where(obs['redshift']==redshift)
		plt.errorbar(obs['lbol'][tt] + loglbol_sun, obs['dphi'][tt],yerr=obs['sig'][tt],fmt='o',color='grey',ecolor='grey',label='Hopkins+ 2007 (Compilation)')
		if z[zz] == 6:
			plt.errorbar([logLbol6B],[logphi6B],xerr=errlogLbol6B,yerr=errlogphi6B,fmt='^',color='k',label='Barger+2003')
	                plt.errorbar([logLbol6F],[logphi6F],xerr=[[logLbol6Fm],[logLbol6Fp]],yerr=[[errlogphi6Fm],[errlogphi6Fp]],fmt='s',color='k',label='Fiore+ 2012')
		if z[zz] == 5:
			l1450 = np.log10(4.4)+mcABconv(obs2['M1450'],c/(0.145e-4))
	                dphi = 10**obs2['logphi']
	                dphip = (2./5.) * (dphi+obs2['sig'])
	                dphim = (2./5.) * (dphi - obs2['sig'])
	                dphi = np.log10((2./5.)*dphi)
	                dphierr = [dphi-np.log10(dphim),np.log10(dphip)-dphi]
	                plt.errorbar(l1450,dphi,yerr=dphierr,fmt='D',color='k',label='McGreer+ 2013')
	if overplot==False:
		plt.title(str(zbinsl[zz[0]])+' < z < '+str(zbinsh[zz[0]]))
     		plt.xlabel(r'log$_{10}$($L_{bol}$ [ergs/s]))',fontsize=30)
     		plt.ylabel(r'log$_{10}$($\phi$ [Mpc$^{-3}$ dex$^{-1}$])',fontsize=30)
        plt.legend(loc='lower left',fontsize=20)
	return

  def run(self):
    dataset = [[float(entry) for entry in data]
               for data in self._get_stripped_file_lines()]
    if not data:
      return None

    # Plot aesthetics.
    fig = plt.figure(1)
    plt.title('%s' % FLAGS.title)
    xlabel = FLAGS.xlabel
    if FLAGS.xlog:
      xlabel = 'log( ' + xlabel + ' )'
    plt.xlabel(xlabel)
    plt.ylabel(FLAGS.ylabel)

    data_plts = []
    for i, data in enumerate(dataset):
      ecdf = distributions.ECDF(data)
      if FLAGS.xmax:
        x = np.linspace(0, float(FLAGS.xmax), num=len(data))
      else:
        x = np.linspace(min(data), max(data), num=len(data))
      y = ecdf(x)
      plt.step(x, y, '.-', label=self.filepaths[i])

    xmin, xmax, ymin, ymax = plt.axis()
    plt.axis((xmin, xmax, 0, 1))
    plt.legend(loc='lower right')
    if FLAGS.xlog:
      plt.xscale('log')
    fig.savefig(FLAGS.plot_name + '.png')

    def hist_alarms(self,alarms,title_str='alarms',save_figure=False,linestyle='-'):
        fontsize=15
        T_min_warn = self.T_min_warn
        T_max_warn = self.T_max_warn
        if len(alarms) > 0:
            alarms = alarms / 1000.0
            alarms = np.sort(alarms)
            T_min_warn /= 1000.0
            T_max_warn /= 1000.0
            plt.figure()
            alarms += 0.0001
            bins=np.logspace(np.log10(min(alarms)),np.log10(max(alarms)),40)
            #bins=linspace(min(alarms),max(alarms),100)
            #        hist(alarms,bins=bins,alpha=1.0,histtype='step',normed=True,log=False,cumulative=-1)
            #
            plt.step(np.concatenate((alarms[::-1], alarms[[0]])), 1.0*np.arange(alarms.size+1)/(alarms.size),linestyle=linestyle,linewidth=1.5)

            plt.gca().set_xscale('log')
            plt.axvline(T_min_warn,color='r',linewidth=0.5)
            #if T_max_warn < np.max(alarms):
            #    plt.axvline(T_max_warn,color='r',linewidth=0.5)
            plt.xlabel('Time to disruption [s]',size=fontsize)
            plt.ylabel('Fraction of detected disruptions',size=fontsize)
            plt.xlim([1e-4,4e1])#max(alarms)*10])
            plt.ylim([0,1])
            plt.grid()
            plt.title(title_str)
            plt.setp(plt.gca().get_yticklabels(),fontsize=fontsize)
            plt.setp(plt.gca().get_xticklabels(),fontsize=fontsize)
            plt.show()
            if save_figure:
                plt.savefig('accum_disruptions.png',dpi=200,bbox_inches='tight')
        else:
            print(title_str + ": No alarms!")

def coc_set_up_plot(bincenters, hist, gaussian, variable, threshold = 0, sub_par = ""):
    '''
    Set up the plotting space for the Climatological Outlier Check

    :param array bincenters: bin centres of histogram
    :param array hist: histogram values
    :param array gaussian: parameters of gaussian fit [m, s, n]
    :param str variable: name of variable for title
    :param int threshold: threshold to plot
    :param str sub_par: sub parameter for axis label
    '''   
    import matplotlib.pyplot as plt
    
    plt.clf()
    plt.axes([0.1,0.15,0.85,0.75])
    plot_hist = np.array([0.01 if h == 0 else h for h in hist])  
    plt.step(bincenters, plot_hist, 'k-', label = 'standardised months', where='mid')

    # plot fitted Gaussian
    plot_gaussian = utils.gaussian(bincenters, gaussian)
    plt.plot(bincenters, plot_gaussian, 'b-', label = 'Gaussian fit')

    # sort the labels etc
    plt.xlabel("%s offset (IQR)" % variable)                    
    plt.ylabel("Frequency (%s)" % sub_par)
    plt.gca().set_yscale('log')
    plt.axvline(-threshold-1,c='r')
    plt.axvline(threshold+1,c='r')
    plt.axvline(-threshold,c='orange')
    plt.axvline(threshold,c='orange')
    plt.ylim(ymin=0.1)
    plt.title("Climatological Gap Check - %s - %s" % (sub_par, variable) )        

    return  # coc_set_up_plot

def plot_concurrency(trans_stats, filename):
    sorted_keys = trans_stats.keys()
    sorted_keys.sort()
    # outfilename = filename[0:filename.rfind(".")]+'_concur.csv'
    # outfile = open(outfilename, "w")
    for key in sorted_keys:
        max_time = 65000
        trans_stats[key].sort()
        step = []
        concur = []
        for i in xrange(0, 65000, 10):
            step.append(i)
            pos = 0
            c = 0
            while (pos < len(trans_stats[key])):
                if (trans_stats[key][pos][0] > i):
                    break
                c = trans_stats[key][pos][1]
                pos += 1
            concur.append(c)
            # line = ("%d,%d") % (i,c)
            # outfile.write(line)
        plt.step(step, concur)
        label='dest_'+key
        # plt.xlim(0, max_time)
        # plt.ylim(0, max_conr+1)
    plt.legend()
    outfile = filename[0:filename.rfind(".")]+'_concur.png'
    plt.savefig(outfile)

def _plot_completeness(comw, start_time, end_time):
    comp = np.column_stack([np.hstack([end_time, comw[:, 0], start_time]),
                            np.hstack([comw[0, 1], comw[:, 1], comw[-1, 1]])])
    for row in comp:
        print(row[0], row[1])
    plt.step(comp[:-1, 0], comp[1:, 1], ls='--',
             where="post", linewidth=3, color='brown')

def plotCDF(list_elements, sim_run, repetition, actual_capacity, TOFILE):
    N = len(list_elements)
    x = numpy.sort(list_elements)
    y = numpy.linspace(0, 1, N)

    plt.step(x, y * 100)
    plt.title('Simulation %d-%d with capacity of %d' %
              (sim_run + 1, repetition + 1, actual_capacity))
    plt.xlabel('average time in system [tu]')
    # displaying tick range. Adjustments needed for long run simulations!
    plt.ylim(0, 100)
    if (max(x) >= 500):
        tickrange = 50
    elif (500 > max(x) >= 100):
        tickrange = 25
    elif (100 > max(x) >= 20):
        tickrange = 5
    else:
        tickrange = 1

    plt.xticks(range(0, int(max(x))+2, tickrange))
    plt.yticks(range(0, 100, 10))
    plt.ylim(ymax=101)
    plt.grid()

    if (TOFILE):
        plt.tight_layout()
        plt.savefig('./figs/cdf/time_in_sys/simulation%d-%d.pdf' %
                    (sim_run + 1, repetition + 1), dpi=(150), format='pdf')
    else:
        plt.show()
    plt.close()

def plotVectors(json, coincidence=0, name="", minTime=30, maxTime=1700,
                labels=False):
    delims = readDelims()

    macAddresses = {}
    json["packets"].apply(lambda row: addPacket(row, macAddresses), axis=1)

    sum = np.zeros(json["last"], dtype=np.int)
    for val in macAddresses.values():
        diff = val[1] - val[0]
        if ((diff > minTime) and (diff < maxTime) and (val[2] >= coincidence)):
            sum[val[0]: val[1]] += 1

    plot.xlabel('Seconds since start')
    plot.ylabel('Devices', color='b')
    plot.step(range(int(json["last"])), sum.tolist())

    actual = [stop['actual'] for stop in delims]
    stopx = [stop['start'] for stop in delims]
    ax2 = plot.twinx()
    ax2.step(stopx, actual, 'r', where="post")
    (t1.set_color for t1 in ax2.get_yticklabels())
    ax2.set_ylabel('Riders', color='r')
    plot.ylim(0, 30)
    plot.xlim(0, json["last"])

    if(labels):
        for stop in delims:
            annotate(stop["code"], stop["start"], stop["actual"], 10, 10)

    makeWidePlot("bus", "vectors")

    plot.show()

	def line_flux(self,EL_wave,wave_range=None,plot=False):
		''' 
		-----PURPOSE-----
		Extract the flux from the emission line from the 'SPEC1D_NORMIVAR' flux density
		-----INPUT-------
		EL_wave 		Cental wavelength of the emission line [Angstroms]
		
		wave_range	    [min,max] where min and max are the wavelenghts in Angstroms
						that start and end the range over which you want to fit. 
						By default (None) will use [EL_wave-25, EL_wave+25]
		
		plot 			if True, will plot the flux density and the 1D Gaussian fit 
		-----OUTPUT------
		line_flux 		In erg/cm^2/s
		'''
		if wave_range == None:
			wave_range = [EL_wave-25,EL_wave+25]
		waves_line, g = self.fit_gaussian(EL_wave=EL_wave,wave_range=wave_range)
		line_flux = sum(g(waves_line))
		# print "Line flux is %.4g erg/s/cm^2" % line_flux
		wave,spec1D_calibrated = self.extract_normalized_spectrum(spectype='SPEC1D_NORMIVAR')
		line_mask = np.logical_and(wave >= min(wave_range),wave <= max(wave_range))
		waves_line = wave[line_mask]
		spec1D_line = spec1D_calibrated[line_mask]
		if plot:
			plt.step(waves_line,spec1D_line,color='r',label='flux density')
			plt.plot(waves_line,g(waves_line),color='b',label='Gaussian fit')
			plt.legend()	
		return line_flux

def plot_hists_lines(feature_list, title_list, fname, colours, buckets, plot_hist,
               main_title_prefix, x1, x2, xlabel, ylabel, obwiednia=False):
    '''
    Plot a number of histograms on the same figure. 
    plot_hist - function responsible for drawing a single histogram, which is passed
        list of values, colour, bins list
        
    '''
    print "[plot_hists_lines]"
    bins = numpy.linspace(min([min(x) for x in feature_list]), 
                          max([max(x) for x in feature_list]), buckets)

    from itertools import izip
    for feats, col in izip(feature_list, colours):
        if obwiednia:
            h1,edges1 = numpy.histogram(feats, bins, normed=False)
            pyplot.step(edges1[1:], h1, color=col)
        else:
            plot_hist(feats, col, bins)
    full_title = []
    if main_title_prefix:
        full_title.append(main_title_prefix)
    #for title, col in izip(title_list, colours):
    #    full_title.append(title+" has colour: "+col+"  ")
    
    pyplot.title(" ".join(full_title))
    #pyplot.axvline(x1, 0, 7000, color="black")
    #pyplot.axvline(x2, 0, 7000, color="black")
    frame1 = pyplot.gca()
    #frame1.get_yaxis().set_visible(False)
    pyplot.xlabel(xlabel)
    pyplot.ylabel(ylabel)
    
    pyplot.savefig(fname)
    pyplot.close()