docs: def prod_fac_massgrid_A_1(self,weighted=True,log=True,runs=[],isotopes=[],elements=[],cycles=[],plot_set_diagram=True,color=['r','b','g','k'],marker_type=['o','p','s','D'],line_style=['--','-','-.',':'],markersize=[6,6,6,6],line_width=[14,14,14,14],title='',withlabel=True,label='',plot_lines=True,exp_only=False,pre_exp=False,delay=True,exp_dir='',fontsizelabel='x-small',iniabupath='/astro/critter/critter/PPN/forum.astro.keele.ac.uk/frames/mppnp/USEEPP/iniab2.0E-02GN93.ppn',xlim=[0,0]):
  		'''
  			Plots behaviour star mass dependent yields of different isotopes - specify dirs at the beginning
  			Beware  of different z,
  
  			runs : If array is empty function uses all available directories.					
  				
  			!! If len of elements longer than zero, than isotopes will be ignored and elements used!!
  
  		'''
  		runs=[]
  
  		HDF5_surf=[]
  		HDF5_out=[]
  		HDF5_restart=[]
  		for i in range(len(self.run_dirs_name)):
  			HDF5_surf.append(self.runs_H5_surf[i])
  			HDF5_out.append(self.runs_H5_out[i])
  			HDF5_restart.append(self.runs_H5_restart[i])
  		legend_k=0
  		sefiles=[]
  		sefiles_hout=[]
  		sefiles_restart=[]
  		for i in range(len(HDF5_surf)):
  				sefiles.append(se(HDF5_surf[i]))		
  				sefiles_hout.append(se(HDF5_out[i]))
  				sefiles_restart.append(se(HDF5_restart[i]))
  
  
  		import utils as u
  		iniabu=u.iniabu(iniabupath)
  		x_iniabu=[]
  
  		
  		isotopes11=sefiles[0].se.isotopes
  		if len(isotopes)>0:
  			if isotopes[0]=='allstable':
  				isotopes=get_stable(isotopes11,get_elements=False)
  				color=len(isotopes)*[color[0]]
  				marker_type=len(isotopes)*[marker_type[0]]
  				line_style=len(isotopes)*[line_style[0]]
  				markersize=len(isotopes)*[markersize[0]]
  				line_width=len(isotopes)*[line_width[0]]
  				print isotopes
  			x_iniabu=iniabu.iso_abundance(isotopes)
  			mass_numbers=iniabu.a
  			charge=iniabu.z
  
  			
  
  
  		if len(elements)>0:
  			if elements[0]=='allstable':
  				elements=get_stable(isotopes11,get_elements=True)
  				color=len(elements)*[color[0]]
  				marker_type=len(elements)*[marker_type[0]]
  				line_style=len(elements)*[line_style[0]]
  				markersize=len(elements)*[markersize[0]]
  				line_width=len(elements)*[line_width[0]]
  				print elements
  			isotopes=get_stable(isotopes11,get_elements=False)
  			import re
  			x_iniabu=[0]*len(elements)
  			x_iniabu_names=[]
  			for k in range(len(iniabu.names)):
  				iso=iniabu.names[k].replace(' ','')
  				charge1=iniabu.z
  				mass_numbers1=iniabu.a
  				match = re.match(r"([a-z]+)([0-9]+)",iso, re.I)
  				ele1=match.groups()[0].upper()
  				iso=ele1+'-'+match.groups()[1]
  				for r in range(len(elements)):
  					if elements[r] == ele1:
  						x_iniabu_names.append(iso)
  						x_iniabu[r]+=iniabu.iso_abundance(iso)
  			print 'test ouptput',x_iniabu_names	
  
  		#import utils as u
  
  		#iniabu=u.iniabu(iniabupath)
  
  		#iniabu.iso_abundance(isotopes)
  
  
  		####dealign with explosion pinput
                  #test to take right dir:
                  #
                  if delay:
                          exp_type='delay'
                  else:
                          exp_type='rapid'
                  slist = os.listdir(exp_dir)
                  expr = re.compile(exp_type)
                  slist=(filter(expr.search,slist))
  
  
  		exp_runs_H5_restart=[]
  		sefiles_exp=[]
  
  
  		#to dinstiungish between pre exp and exp sources
  		if pre_exp==False:
  
  			slist = os.listdir(exp_dir)
  			expr = re.compile(exp_type)
  			slist=(filter(expr.search,slist))
  
  			for element in slist:
  				run_path=exp_dir+'/'+element
  				if not os.path.isdir(run_path):
  					continue
  				if os.path.isdir(run_path+"/H5_restart"):
  					sefiles1 = os.listdir(run_path+"/H5_restart")
  					if (filter(expr.search,sefiles1)) <1:
  						print "Warning: No hdf5 restart files found in "+run_path+"/H5_restart"
  					else:
  						exp_runs_H5_restart.append(run_path+"/H5_restart")
  						sefiles_exp.append(se(run_path+"/H5_restart"))
  		else:
  	
  			exp_runs_H5_restart=HDF5_restart
  			for k in range(len(HDF5_restart)):
  				sefiles_exp.append(se(HDF5_restart[k]))
  
  
  		z_index_files=[]
  		z_values=[]
  		j=-1
  		for i in range(len(HDF5_surf)):
  			j+=1
  			star_z=sefiles[i].get("zini")
  			if star_z not in z_values:
  				z_values.append(star_z)
  				z_index_files.append([])
  			z_index_files[ z_values.index(star_z)].append(i)
  		max_yield=[]
  		color_iso=-1
  		yields_1=[]
  		t=0
  		if len(elements)>0:
  			isotopes=elements
  		legend_k=0
  		for w in range(len(z_index_files)):
  			star_mass_array=[]
  			yields=[]
  			legend_k+=1
  			#iso_yields=[]
  			#iniabu_yields_folded=[]
  			production_factor=[]
  			for i in range(len(isotopes)):
  				#iso_yields.append(np.zeros(len( z_index_files[w]   )))
  				#iniabu_yields_folded.append(np.zeros(len( z_index_files[w]   )))
  				if exp_only==False:	
  					production_factor.append(np.zeros(len( z_index_files[w]   )))
  				else:
  					production_factor.append([])
  			ttt=0
  			for k in z_index_files[w]:
  				if type(sefiles[k].get("mini")) == np.ndarray:
  					star_mass=sefiles[k].get("mini")[0]
  					star_z=sefiles[k].get("zini")[0]
  				else:
                                          star_mass=sefiles[k].get("mini")
                                          star_z=sefiles[k].get("zini")
  					
  					
  				#star_mass=sefiles[k].get("mini")[0]
  				#star_z=sefiles[k].get("zini")[0]				
  				if cycles[k][1]==-1:
  					endcycle=int(sefiles[k].se.cycles[-1]) #+ cycles[k][2]  #1000
  					endcycle= int( round( endcycle,-3) -2000 )
  				else:
  					endcycle=cycles[k][1]
  				if exp_only==False:
  					star_mass_array.append(star_mass)
  					prod_factor,isotopes_prod_fac,yields,iniabu,mass_frac_ini,remn_mass =self.weighted_yields(sefiles[k],sefiles_hout[k],sefiles_restart[k],isotopes,elements,cycles[k][0],endcycle,cycles[k][2])
  					#print 'yield output###################################:'
  					#print 'wind: ',yields,prod_factor,iniabu
  					for tt in range(len(prod_factor)):
                                          	prod_factor[tt]=yields[tt]/(star_mass-remn_mass)    #iniabu[tt]
  
  
  				else:
  					prod_factor=[]
  					if star_mass<=8:
  						continue
  					for pp in range(len(isotopes)): 
  						production_factor[pp].append([0])
  							
  					star_mass_array.append(star_mass)
  					
  				mass=star_mass
  				metallicity=star_z
  				if mass>8:
  
  					#in case no ppn_exp run dir is available, skip explosion contribution
  					#this is because mass cut is then larger then actual mass of star...no mass lost
  
  					for t in range(len(exp_runs_H5_restart)):
  						if type(sefiles_exp[t].get("mini")) == np.ndarray:
  							mass1=sefiles_exp[t].get("mini")[0]
  							metallicity1=sefiles_exp[t].get("zini")[0]
  						else:
  							mass1=sefiles_exp[t].get("mini")
  							metallicity1=sefiles_exp[t].get("zini")
  
  						###for currupted files
  
  						corrupt=False
  						if type(mass1) == list:
  							corrupt=True
  						if type(metallicity1) ==list:
  							corrupt=True
  						if ((float(mass1) < 1) or  (float(mass1) >70 )) or type(mass1) == type(None):
  							corrupt=True
  						if ((float(metallicity1) < 0.00001) or  (float(metallicity1) >1. )) or type(metallicity1) == type(None):
  							corrupt=True
  
  
  						if corrupt == True:
  							metallicity1=0.02
  							rundir11=exp_runs_H5_restart[t].split('/')[-2]
  							mass1=float(rundir11.split('.')[0][1:])
  							#print 'corruption, assinge new values',mass1,metallicity1
  						if mass1 == mass and metallicity1 == metallicity:
  							#sefiles_re2=sefiles_exp[t]
  							#calculate exp part
  							import nugridse as mp
  							reload(mp)
  							#Have to find the right cycle in restart file
  							#.se.cycles does not work when reading all files	
  							#therefore identify file with last cycle and read it:
  							refiles=sefiles_exp[t].se.files
  							tosort=[]
  							refiles1=[]
  							for i in range(len(refiles)):
  								cycle_file=refiles[i][-18:-11]
  								if cycle_file.isdigit():
  									tosort.append(int(cycle_file))
  									refiles1.append(refiles[i])
  							idx_sorted=sorted(range(len(tosort)), key=lambda k: tosort[k])
  							idx_file = idx_sorted[-1]
  							sefiles_re_cycle=mp.se(exp_runs_H5_restart[t],refiles1[idx_file])
  							if len(sefiles_re_cycle.se.cycles)==0:
  								#in the case there is not cycle in the last file...set1.2 M2-
  								sefiles_re_cycle=mp.se(exp_runs_H5_restart[t],refiles1[idx_sorted[-2]])
  							print 'tewststest',exp_cycle_inp
  							if exp_cycle_inp>0:
  								sefiles_re_cycle=mp.se(exp_runs_H5_restart[t],exp_cycle_inp)
  							print 'Use custom exp cycle ', exp_cycle_inp
  							#print 'identifies as ',exp_runs_H5_restart[t],'with file',refiles1[idx_file]
  							isotopes_exp,yields_exp,iniabu_exp,mass_cut_exp,first_cycle_exp=self.weighted_yields_explosion(mp,sefiles_re_cycle,sefiles[k],isotopes=isotopes,elements=elements,delay=delay)
  
  							if exp_only==False:	
  								#add pre-exp yields
  								for tt in range(len(prod_factor)):
  									#print 'star above 8M'
  									#print isotopes_prod_fac[tt],yields[tt],yields_exp[tt]
  									#the iniabu from the explosion is not added
  									#because it is already in the iniabu from wind
  									prod_factor[tt]=(yields[tt]+yields_exp[tt])/(star_mass-mass_cut_exp)     #(iniabu[tt])
  									#print 'yields wind:',yields[tt]
  									#print 'prodfac wind:',(yields[tt]/iniabu[tt])
  									#print 'wind+exp yields : ',yields[tt]+yields_exp[tt]	
  							else:
  								for tt in range(len(isotopes_exp)):		
  									prod_factor.append((yields_exp[tt])/( (star_mass-mass_cut_exp) - (star_mass-remn_mass) ))   #(iniabu_exp[tt]))
  							remn_mass=mass_cut_exp
  							#print 'exp yield',yields_exp
  							#print 'prodfac exp:',(yields_exp[tt]/iniabu_exp[tt])
  							#print 'wind+exp prodfac:',prod_factor[tt]
  							break		
  
  
  		
  
  				#print 'iso exp',len(isotopes_exp)
  				#prod_factor=np.log10(prod_factor)
  				print len(isotopes),len(prod_factor),len(production_factor)
  				#Normalization to solar abundance
  				for i in range(len(isotopes)):
                                          production_factor[i][ttt]=prod_factor[i]/x_iniabu[i]
  				ttt+=1
  			
  			###plotting
  			#if plot_set_diagram==True:
  			#if plot_set_diagram==True:
  			#	fig_2=plt.figure(isotopes[h])
  			#	#fig_2=plt.figure()
  			#	ax_2 = fig_2.add_subplot(1,1,1)
  			#	plt.rcParams.update({'font.size': 16})
  			#	plt.rc('xtick', labelsize=16)
  			#	plt.rc('ytick', labelsize=16)
  			#	if log==True:
  			#		ax_2.set_yscale('log')
                          color_iso=0
                          legend_k=0
  
  
  			for h in range(len(isotopes)):
  				#yield_1=iso_yields[i]
  				#mass_1=star_mass_array
  				yield_1=[]
  				mass_1=[]			
  				iniabu_folded=[]
  				prod_fac_sorted=[]	
                 			indices_array=sorted(range(len(star_mass_array)),key=lambda x:star_mass_array[x])
                  		for i in indices_array:
  				#	yield_1.append(iso_yields[h][i])
  					mass_1.append(star_mass_array[i])
  				#	iniabu_folded.append(iniabu_yields_folded[h][i])
  					prod_fac_sorted.append(production_factor[h][i])
  				####Plotting prodfactor
  				#plt.figure(fig_2.number)
  
  				fig_2=plt.figure(isotopes[h])
  				ax_2 = fig_2.add_subplot(1,1,1)
  				plt.rcParams.update({'font.size': 16})
  				plt.rc('xtick', labelsize=16)
  				plt.rc('ytick', labelsize=16)
  				if log==True:
  					ax_2.set_yscale('log')
  				ax_2.legend()
  				ax_2.set_xlabel("M/M$_{\odot}$",fontsize=16)
  				ax_2.minorticks_on()
  				ax_2.set_ylabel("Overproduction factor",fontsize=16)
  				ax_2.set_title(title)
  	
  				if len(label)>0:
  					label=", "+label
  				if len(elements)==0:
  					plot_quantity=isotopes[h]
  					plot_quantity="$^{"+plot_quantity.split("-")[1]+"}$"+plot_quantity.split("-")[0]
  				else:
  					plot_quantity=elements[h]
  				if withlabel==True:
                                  	plt.plot(mass_1,prod_fac_sorted,marker=marker_type[legend_k],color=color[color_iso],markersize=markersize[legend_k],mfc=color[color_iso],linewidth=line_width[legend_k],linestyle=line_style[color_iso],label=plot_quantity+" , Z="+str(star_z)+label  )				
  				else:
                                  	plt.plot(mass_1,prod_fac_sorted,marker=marker_type[legend_k],color=color[color_iso],markersize=markersize[legend_k],mfc=color[color_iso],linewidth=line_width[legend_k],linestyle=line_style[color_iso])				
  				
  
  				m_max=max(star_mass_array)+2
  				
  				if plot_lines==True:
  					plt.plot([0,m_max],[1,1],"k--",linewidth=3)
  					plt.plot([0,m_max],[2,2],"k--",linewidth=1)
  					plt.plot([0,m_max],[0.5,0.5],"k--",linewidth=1)				
  	
  				color_iso+=1
  				legend_k+=1
  		#####
  		x_imf=[]
  		y_imf=[]
  		m_max=max(star_mass_array)+2.