new models and update to handle MaNGA DR17

python/firefly_models.py

@@ -114,7 +114,7 @@ def __init__(self, specObs, outputFile, cosmo, models = 'MaStar', model_libs = [
 		self.use_downgraded_models = use_downgraded_models
 		self.write_results = write_results
 		self.flux_units = flux_units
-		if self.models == 'm11':
+		if (self.models == 'm11') or (self.models == 'm11-sg'):
 			for m in self.model_libs:
 				if m == 'MILES':
 					self.deltal_libs.append(2.55)
@@ -127,7 +127,7 @@ def __init__(self, specObs, outputFile, cosmo, models = 'MaStar', model_libs = [
 
 		elif self.models =='MaStar':
 			model_path = os.environ['STELLARPOPMODELS_DIR']
-			ver='v0.2'
+			ver='v1.1'
 			hdul=pyfits.open(model_path+'/MaStar_SSP_'+ver+'.fits.gz')
 			r_model=hdul[2].data[1,:]
 			# This provides R=lamba/delta_lambda as numpy ndarray. The params deltal_libs and deltal should probably be renamed. 
@@ -181,6 +181,119 @@ def get_model(self, model_used, imf_used, deltal, vdisp, wave_instrument, r_inst
 			and returns it as well
 
 		"""
+
+		if self.models == 'm11-sg':
+			first_file  = True
+			model_files = []
+			#print('yes we are in here')
+			#stop
+#			if self.use_downgraded_models :
+#				if model_used == 'MILES_UVextended' or model_used == 'MILES_revisedIRslope':
+#					model_path 		= join(os.environ['STELLARPOPMODELS_DIR'],'SSP_M11_MILES_downgraded','ssp_M11_' + model_used+ '.' + imf_used)
+#				else:
+#					model_path 		= join(os.environ['STELLARPOPMODELS_DIR'],'SSP_M11_'+ model_used + '_downgraded', 'ssp_M11_' +model_used +'.' + imf_used)
+#			else:
+#				if model_used == 'MILES_UVextended' or model_used == 'MILES_revisedIRslope':
+#					model_path 		= join(os.environ['STELLARPOPMODELS_DIR'],'SSP_M11_MILES', 'ssp_M11_'+model_used+'.'+imf_used)
+#				else:
+			model_path 		= join(os.environ['STELLARPOPMODELS_DIR'],'SSP_M11_'+model_used+'_SG' ,'ssp_M11_' +model_used +'.' + imf_used)
+
+
+			# Constructs the metallicity array of models :
+			all_metal_files = sorted(glob.glob(model_path+'*'))
+			#print(model_path)
+			#print(all_metal_files)
+			#stop	
+			## # print all_metal_files
+			metal_files 	= []
+			metal 	    = [] #[-2.25, -1.35, -0.33, 0, 0.35]
+			for z in range(len(all_metal_files)):
+				zchar = all_metal_files[z][len(model_path):]
+				if zchar == 'z001.sg':
+					znum = 10**(-0.33) 
+				elif zchar == 'z002.sg':
+					znum = 10**(0) 
+				elif zchar == 'z004.sg':
+					znum = 10**(0.35)
+				elif zchar == 'z0001.bhb.sg':
+					#znum = -1.301
+					znum = 10**(-1.35) #10**-1.301
+				elif zchar == 'z0001.rhb.sg':
+					#znum = -1.302
+					znum = 10**(-1.35) #10**-1.302
+				#elif zchar == 'z10m4.bhb':
+					#znum = -2.301
+					#znum = 10**(-2.25) #10**-2.301
+				#elif zchar == 'z10m4.rhb':
+					#znum = -2.302
+					#znum = 10**(-2.25) #10**-2.302
+				#elif zchar == 'z10m4':
+					#znum = -2.300
+					#znum = 10**(-2.25) #10**-2.300
+				elif zchar == 'z0p25.sg':
+					znum = 10**0.25
+				elif zchar == 'zm0p7.bhb.sg':
+					znum = 10**-0.7
+				elif zchar == 'zm0p7.rhb.sg':
+					znum = 10**-0.7
+				elif zchar == 'zm1p0.bhb.sg':
+					znum = 10**-1.0
+				elif zchar == 'zm1p0.rhb.sg':
+					znum = 10**-1.0
+				else:
+					raise NameError('Unrecognised metallicity! Check model file names.')
+
+				if znum>10**(self.Z_limits[0]) and znum<10**(self.Z_limits[1]):
+					metal_files.append(all_metal_files[z])
+					metal.append(znum)
+			#print(metal_files)
+			#stop
+			# constructs the model array
+			model_flux, age_model, metal_model = [],[],[]
+			for zi,z in enumerate(metal_files):
+				# print "Retrieving and downgrading models for "+z
+				model_table = pd.read_table(z,converters={'Age':np.float64}, header=None ,usecols=[0,2,3], names=['Age','wavelength_model','flux_model'], delim_whitespace=True)
+				age_data = np.unique(model_table['Age'].values.ravel())
+#				print(age_data)
+#				stop
+				for a in age_data:
+					logyrs_a = trylog10(a)+9.0
+					## print "age model selection:", self.age_limits[0], logyrs_a, self.age_limits[1]
+					if (((10**(logyrs_a-9)) < self.age_limits[0]) or ((10**(logyrs_a-9)) > self.age_limits[1])):
+						continue
+					else:
+						spectrum = model_table.loc[model_table.Age == a, ['wavelength_model', 'flux_model'] ].values
+						wavelength_int,flux = spectrum[:,0],spectrum[:,1]
+
+						# converts to air wavelength
+						if self.data_wave_medium == 'vacuum':
+							wavelength = airtovac(wavelength_int)
+						else:
+							wavelength = wavelength_int
+
+						# downgrades the model
+						if self.downgrade_models:
+							mf = downgrade(wavelength,flux,deltal,self.specObs.vdisp, wave_instrument, r_instrument)
+						else:
+							mf = copy.copy(flux)
+
+						# Reddens the models
+						if ebv_mw != 0:
+							attenuations = unred(wavelength,ebv=0.0-ebv_mw)
+							model_flux.append(mf*attenuations)
+						else:
+							model_flux.append(mf)
+
+						age_model.append(a)
+						metal_model.append(metal[zi])
+						first_model = False
+			#print(wavelength)
+			#stop
+
+			# print "Retrieved all models!"
+			self.model_wavelength, self.model_flux, self.age_model, self.metal_model = wavelength, model_flux, age_model, metal_model
+			return wavelength, model_flux, age_model, metal_model
+
 		# first the m11 case
 		if self.models == 'm11':
 			first_file  = True
@@ -276,7 +389,7 @@ def get_model(self, model_used, imf_used, deltal, vdisp, wave_instrument, r_inst
 
 						# downgrades the model
 						if self.downgrade_models:
-							mf = downgrade(wavelength,flux,deltal,self.vdisp_round, wave_instrument, r_instrument)
+							mf = downgrade(wavelength,flux,deltal,self.specObs.vdisp, wave_instrument, r_instrument)
 						else:
 							mf = copy.copy(flux)
 
@@ -300,7 +413,7 @@ def get_model(self, model_used, imf_used, deltal, vdisp, wave_instrument, r_inst
 		elif self.models =='MaStar':
 
 			model_path = os.environ['STELLARPOPMODELS_DIR']
-			ver = 'v0.2'
+			ver = 'v1.1'
 
 			lib = model_used
 			if imf_used == 'kr':
@@ -323,10 +436,8 @@ def get_model(self, model_used, imf_used, deltal, vdisp, wave_instrument, r_inst
 
 			wavelength=hdul[2].data[0,:]
 
-			if (lib=='Th-MaStar'):
+			if (lib=='gold'):
 				fluxgrid=hdul[3].data
-			if (lib=='E-MaStar'):
-				fluxgrid=hdul[4].data
 
 			sidx = np.where(s==slope)[0][0]
 
@@ -345,7 +456,7 @@ def get_model(self, model_used, imf_used, deltal, vdisp, wave_instrument, r_inst
 
 					# downgrades the model
 					if self.downgrade_models:
-						mf = downgrade(wavelength,flux,deltal,self.vdisp_round, wave_instrument, r_instrument)
+						mf = downgrade(wavelength,flux,deltal,self.specObs.vdisp, wave_instrument, r_instrument)
 					else:
 						mf = copy.copy(flux)
 

python/firefly_setup.py

@@ -157,7 +157,7 @@ def openSingleSpectrum(self, wavelength, flux, error, redshift, ra, dec, vdisp,
 		#print(self.ebv_mw)	
 
 
-	def openMANGASpectrum(self, path_to_logcube, path_to_dapall, bin_number, plate_number, ifu_number, emlines):
+	def openMANGASpectrum(self, path_to_logcube, path_to_dapall, bin_number, plate_number, ifu_number, emlines,mpl='mpl-9'):
 		"""Loads an observed MaNGA spectrum in.
 		:param path_to_logcube: Must specify the path to logcube (if using MPL5 or higher). Set to 0 otherwise.		
 		"""
@@ -170,8 +170,15 @@ def openMANGASpectrum(self, path_to_logcube, path_to_dapall, bin_number, plate_n
 
 		# Get S/N, right ascension and declination.
 		signal, ra, dec = maps_header['BIN_SNR'].data[x_position,y_position], maps_header[0].header['OBJRA'],maps_header[0].header['OBJDEC']
-		velocity_dispersion = maps_header['STELLAR_SIGMA'].data 		# DO NOT USE VELOCITY DISPERSION CORRECTION!
-		vdisp = velocity_dispersion[x_position,y_position]
+		velocity_dispersion = maps_header['STELLAR_SIGMA'].data 		# DO NOT USE VELOCITY DISPERSION CORRECTION!		
+		velocity_dispersion_correction = maps_header['STELLAR_SIGMACORR'].data[0,:,:]
+
+		if velocity_dispersion[x_position,y_position] > velocity_dispersion_correction[x_position,y_position]:
+			correction = np.sqrt((velocity_dispersion[x_position,y_position])**2-(velocity_dispersion_correction[x_position,y_position])**2)
+			vdisp = correction
+		else:
+			vdisp = 0
+
 
 		# Open LOGCUBE to get the flux, wavelength, and error
 		header = pyfits.open(path_to_logcube)
@@ -182,19 +189,30 @@ def openMANGASpectrum(self, path_to_logcube, path_to_dapall, bin_number, plate_n
 		output_flux = correct_flux - correct_flux_emline
 		correct_inverse_variance = inverse_variance[:, x_position, y_position]
 
+		LSF = header['LSF'].data[:,x_position,y_position]		# LSF given as sigma of Gaussian in Angstrom
+		sig2fwhm        = 2.0 * np.sqrt(2.0 * np.log(2.0))
+		LSF_FWHM = LSF*sig2fwhm
+		RES = wavelength/LSF_FWHM
+
+		self.r_instrument = RES
 		self.error = np.sqrt(1.0/(correct_inverse_variance))
 		self.bad_flags = np.ones(len(output_flux))
 		self.flux = output_flux
 		self.vdisp = vdisp
+
+		if (mpl=='mpl-10') or (mpl=='mpl-11'):
+			ext=2
+		else:
+			ext=1
 
 		dap_all = pyfits.open(path_to_dapall)
-		get = np.where(dap_all[1].data['PLATEIFU']==str(plate_number)+'-'+str(ifu_number))
+		get = np.where(dap_all[ext].data['PLATEIFU']==str(plate_number)+'-'+str(ifu_number))
 		c = const.c.value/1000
 		# Use redshift as measured from the stellar kinematics by the DAP.
-		redshift = dap_all[1].data['STELLAR_Z'][get][0]
+		redshift = dap_all[ext].data['STELLAR_Z'][get][0]
 		# If redshift measurement failed, use redshift estimate from NSA or ancillary programs.
 		if redshift<0:
-			redshift = dap_all[1].data['Z'][get][0]
+			redshift = dap_all[ext].data['Z'][get][0]
 
 		sys_vel = maps_header[0].header['SCINPVEL']
 		bin_vel = maps_header['STELLAR_VEL'].data[x_position,y_position]	
@@ -218,8 +236,9 @@ def openMANGASpectrum(self, path_to_logcube, path_to_dapall, bin_number, plate_n
 		self.bad_flags = self.bad_flags[(self.final_mask==False)]
 
 		# Get Trust flag, object_id, xpos, ypos and instrumental resolution.
-		self.trust_flag, self.objid, self.r_instrument = True, 0, np.loadtxt(os.path.join(os.environ['FF_DIR'],'data/MaNGA_spectral_resolution.txt'))
-		self.r_instrument = self.r_instrument[0:self.r_instrument.shape[0]//2]
+# 		self.trust_flag, self.objid, self.r_instrument = True, 0, np.loadtxt(os.path.join(os.environ['FF_DIR'],'data/MaNGA_spectral_resolution.txt'))
+		self.trust_flag, self.objid= True, 0
+# 		self.r_instrument = self.r_instrument[0:self.r_instrument.shape[0]//2]
 		self.r_instrument = self.r_instrument[(self.final_mask==False)]
 		self.xpos, self.ypos = ra, dec
 

run/firefly_manga.py

@@ -36,18 +36,23 @@
 logcube_dir = 'example_data/manga/'
 maps_dir = 'example_data/manga/'
 output_dir = 'output/manga/'
-dap_file = 'example_data/manga/dapall-v2_4_3-2.2.1.fits'
+dap_file = 'example_data/manga/dapall-v3_1_1-3.1.0.fits'	# For DR15 use dapall-v2_4_3-2.2.1.fits
 suffix = ''
 
 # define plate number, IFU number, bin number
 plate = 8080
 ifu = 12701
 bin_number = 0 #'all' if entire IFU, otherwise bin number
 
+# define MaNGA Product Launch (MPL) to use
+# mpl-11 = DR17 final data realease (recommended)
+# mpl-7 = DR15
+mpl = 'mpl-11'
+
 # masking emission lines
 # defines size of mask in pixels
 # set to value>0 for masking (20 recommended), otherwise 0
-N_angstrom_masked=20 
+N_angstrom_masked=0 
 # set emission lines to be masked, comment-out lines that should not be masked
 emlines = [
 						'He-II',# 'He-II:  3202.15A, 4685.74'
@@ -70,13 +75,13 @@
 					    'Ar-III',#'Ar-III: 7135.67'
 						]
 
-# choose model: 'm11', 'MaStar')
+# choose model: 'm11','m11-sg','MaStar')
 model_key = 'MaStar'
 
 # model flavour
 # m11: 'MILES', 'STELIB', 'ELODIE', 'MARCS'
-# MaStar: 'Th-MaStar', 'E-MaStar'
-model_lib = ['E-MaStar']
+# MaStar: 'gold'
+model_lib = ['gold']
 
 # choose IMF: 'kr' (Kroupa), 'ss' (Salpeter)
 imfs = ['kr']
@@ -113,11 +118,11 @@
 def f(i):
     galaxy_bin_number  = i
     print('Fitting bin number {}'.format(i))
-    output_file = direc+'/spFly-'+str(plate)+'-'+str(ifu)+'-bin'+str(i)+'.fits'
+    output_file = direc+'/spFly-'+str(plate)+'-'+str(ifu)+'-bin'+str(i)
     spec = fs.firefly_setup(maps, milky_way_reddening=milky_way_reddening, \
                                   N_angstrom_masked=N_angstrom_masked,\
                                   hpf_mode=hpf_mode)
-    spec.openMANGASpectrum(logcube, dap_file, galaxy_bin_number, plate, ifu, emlines)
+    spec.openMANGASpectrum(logcube, dap_file, galaxy_bin_number, plate, ifu, emlines, mpl=mpl)
 
     age_min = age_limits[0]
     if type(age_limits[1])==str:
@@ -156,8 +161,9 @@ def f_mp(unique_bin_number):       #multiprocessing function
 print('Plate = {}, IFU = {}'.format(plate,ifu))
 
 # Set MAPS and LOGCUBE paths.
-logcube = os.path.join(logcube_dir,'manga-'+str(plate)+'-'+str(ifu)+'-LOGCUBE-VOR10-GAU-MILESHC.fits.gz')
-maps = os.path.join(maps_dir,'manga-'+str(plate)+'-'+str(ifu)+'-MAPS-VOR10-GAU-MILESHC.fits.gz')
+# For DR17 use defaul, for DR15 use '[...]VOR10-GAU-MILESHC.fits.gz'
+logcube = os.path.join(logcube_dir,str(plate),str(ifu),'manga-'+str(plate)+'-'+str(ifu)+'-LOGCUBE-VOR10-MILESHC-MASTARSSP.fits.gz')
+maps = os.path.join(maps_dir,str(plate),str(ifu),'manga-'+str(plate)+'-'+str(ifu)+'-MAPS-VOR10-MILESHC-MASTARSSP.fits.gz')
 
 # Create output path if it doesn't exist.
 direc = os.path.join(output_dir,str(plate),str(ifu))
@@ -170,8 +176,9 @@ def f_mp(unique_bin_number):       #multiprocessing function
 print('Number of bins = {}'.format(len(unique_bin_number)))
 
 if bin_number=='all':
-  N = f_mp(unique_bin_number)
+    N = f_mp(unique_bin_number)
+
 else:
-  f(bin_number)
+    f(bin_number)
 
 print('Time to complete: ', (time.time())-t0)