Added some improvements to the Python scripts for handling spectropolarimetric results (31460c79) · Commits · polarization / boundary_layer

build_model/make_model.py

+167 −108

Original line number	Diff line number	Diff line
		@@ -67,6 +67,9 @@ for inputval in sys.argv:
		if inputval.startswith('emin_model') == True:
		emin_model = float(inputval[inputval.rfind('=') + 1:])

		if inputval.startswith('emax_model') == True:
		emax_model = float(inputval[inputval.rfind('=') + 1:])

		if inputval.startswith('rebin') == True:
		rebin = str(inputval[inputval.rfind('=') + 1:])

		@@ -76,6 +79,9 @@ for inputval in sys.argv:
		if inputval.startswith('pdf') == True:
		pdf = str(inputval[inputval.rfind('=') + 1:])

		if inputval.startswith('plot') == True:
		plot = str(inputval[inputval.rfind('=') + 1:])

		try:
		runtable
		except:
		@@ -121,7 +127,7 @@ except:
		try:
		ymin
		except:
		ymin=-5
		ymin=None

		try:
		ymax
		@@ -134,22 +140,20 @@ except:
		timeon='n'

		try:
		rebin

		rebin = rebin.lower()
		if rebin == 'y':
		rebin = True
		elif rebin == 'n':
		rebin = False
		else:
		print('rebin parameter not set correctly')
		exit()
		raise ValueError("Parameter 'rebin' must be 'y' or 'n'")
		except:
		rebin=False
		rebin = True

		try:
		factor
		except:
		factor=4
		factor=8

		try:
		pdf
		@@ -162,6 +166,10 @@ try:
		except:
		pdf=False

		try:
		plot
		except:
		plot='y'

		#=====================================================================================
		#Energy bins
		@@ -170,18 +178,28 @@ except:
		#====================================================================================

		energy_bins=257
		emin=0.3
		try:
		emin
		except:
		emin=0.1

		try:
		emax
		except:
		emax=70
		emax=150

		try:
		emin_model
		except:
		emin_model=emin

		try:
		emax_model
		except:
		emax_model=emax



		#=============================================================================

		logbins = np.logspace(np.log10(emin), np.log10(emax), num=energy_bins)
		@@ -190,7 +208,7 @@ ene_center=np.asarray([round((logbins[ii]+logbins[ii+1])/2,4) for ii in range(en
		delta_energy=np.asarray([round(logbins[ii+1]-logbins[ii],4) for ii in range(energy_bins-1)])


		logbins_model = np.logspace(np.log10(emin_model), np.log10(emax), num=energy_bins)
		logbins_model = np.logspace(np.log10(emin_model), np.log10(emax_model), num=energy_bins)


		#============================================================================
		@@ -231,6 +249,7 @@ u_centers = 0.5 * (u_bin_edges[:-1] + u_bin_edges[1:])
		#u_centers = u_centers[::-1]



		#for ii in range(len(u_bin_edges)-1):
		# print('Bin %d [%4.2f : %4.2f]' % (ii, u_bin_edges[ii], u_bin_edges[ii+1]))

		@@ -247,7 +266,7 @@ def simul2fits(ene_center=None, delta_ene=None, counts=None, err_counts=None, mc
		emin=ene_center[ii]-delta_ene[ii]/2
		emax= ene_center[ii] + delta_ene[ii] / 2

		f.write('%4.3f %4.3f %4.3f %4.3f\n' % (emin, emax, counts[ii], err_counts[ii]))
		f.write('%5.4f %5.4f %4.3f %4.3f\n' % (emin, emax, counts[ii], err_counts[ii]))

		f.close()

		@@ -284,52 +303,87 @@ def table_init(init_val=None, array_values=None, name=None, table=None):

		#================================================================================================
		def write_table_value(latmax=None, ktbb=None, kte=None, tau=None, u_obs=None, spin=None, mc_file=None,
		xspec=None, fitsfile=None, paramvalues=None, emin_model=None, rebin=None, factor=None):
		xspec=None, fitsfile=None, paramvalues=None,
		emin_model=None, emax_model=None, rebin=None, factor=None, energy_center=None, delta_energy=None):

		print('==============================================')
		print('Reading file: %s' % (mc_file))

		C_matrix, Q_matrix = read_mc_simulations(mc_file=mc_file, n_i=20, n_j=256)
		print('Reading file:\n%s' % (mc_file))
		print('==============================================')


		#Units of C_filtered are counts/keV because the fit is performed by dividing counts
		#of the simulation for the corresponding delta_energy bin
		C_model_kev = smooth_matrix_spectrum(specpol_array=C_matrix, energy_center=ene_center,
		delta_energy=delta_energy, ktbb=ktbb, kte=kte, tau=tau)
		#print(Q_matrix.shape, len(ene_center), len(delta_energy))

		print('Parameters value', ktbb, kte, tau)

		#Change of sign because of convenction: for MC, negative Q values corresponds to PA parallel to the slab
		#while the adopted convection here is the opposite

		Q_matrix=-Q_matrix
		matrix_data=MatrixData(mcfile=mcfile, energy_center=energy_center, delta_energy=delta_energy)
		matrix_data.read_values()
		matrix_data.rebin_data(factor=factor)

		# ===============================================================================
		# Write the I,Q,U spectra
		# ===============================================================================
		C_model_kev = smooth_matrix_spectrum(specpol_array=matrix_data.C_matrix, energy_center=energy_center,
		delta_energy=delta_energy, ktbb=ktbb, kte=kte, tau=tau, emax=150)

		PD_filtered = smooth_matrix_pd(Q_matrix=Q_matrix, C_matrix=C_matrix,
		energy_center=ene_center, delta_energy=delta_energy,
		emin=emin, emax=30, rebin=rebin, factor=factor)
		matrix_data.C_model_kev = C_model_kev

		'''
		Q_model=smooth_matrix_Q_poly(Q_matrix=matrix_data.Q_matrix_rebin, C_matrix=matrix_data.C_matrix_rebin,
		energy_center=energy_center,
		delta_energy=delta_energy,
		energy_center_rebin=matrix_data.energy_center_rebin,
		delta_energy_rebin=matrix_data.delta_energy_rebin,
		emin=emin, emax=emax)

		matrix_data.Q_model = Q_model

		PD_filtered=matrix_data.Q_model/matrix_data.C_model_kev
		matrix_data.Q_model = PD_filtered




		plt.xlim(0.1, 30)
		plt.xscale('log')
		plt.yscale('log')
		plt.plot(energy_center[energy_mask], matrix_data.Q_model[5][energy_mask]*100, color='blue')
		plt.plot(energy_center[energy_mask], matrix_data.C_model_kev[5][energy_mask] * 100, color='orange')
		plt.show()
		'''

		#===================================================
		#Routing performing fit of the Q/C data
		#===================================================

		PD_filtered = smooth_matrix_pd(Q_matrix=matrix_data.Q_matrix_rebin, C_matrix=matrix_data.C_matrix_rebin,
		energy_center=energy_center,
		delta_energy=delta_energy,
		energy_center_rebin=matrix_data.energy_center_rebin,
		delta_energy_rebin=matrix_data.delta_energy_rebin,
		emin=emin, emax=35)

		matrix_data.Q_model = PD_filtered

		#print(PD_filtered.shape)
		#print(C_matrix.shape, Q_matrix.shape, C_model_kev.shape, PD_filtered.shape)

		if C_matrix.shape[0] != len(u_centers):
		if matrix_data.C_matrix_rebin.shape[0] != len(u_centers):
		print('Error: number of u angles from the C_matrix is %d' % (C_matrix.shape[0]))
		print('while the array u_center has %d values ' % (len(u_centers)))
		exit()

		Norm_flux = normalize_counts(counts_array=C_matrix, energy_center=ene_center, delta_energy=delta_energy,

		Norm_flux = normalize_counts(counts_array=matrix_data.C_matrix, energy_center=ene_center, delta_energy=delta_energy,
		u_array=u_centers, delta_u=delta_u, ktbb=float(ktbb))

		I_matrix = intensity_matrix(counts_array=C_model_kev, energy_center=ene_center, delta_energy=delta_energy,
		delta_u=delta_u, norm_value=Norm_flux)


		if fitsfile is None:
		fitsfile = 'intensity_pd.fits'

		try:
		create_fits_image(I_matrix=I_matrix, PD_matrix=PD_filtered, emin=emin, emax=emax, filename=fitsfile)
		create_fits_image(I_matrix=I_matrix, PD_matrix=PD_filtered, energy_center=energy_center, filename=fitsfile)
		except:
		print('Error while creating the fits file %s ' % (fitsfile))
		os.sys.exit(1)
		@@ -349,7 +403,8 @@ def write_table_value(latmax=None, ktbb=None, kte=None, tau=None, u_obs=None, sp
		(float(latmax), float(ktbb), float(kte), float(tau), iobs, spin))


		call_model(fitsfile=fitsfile, iobs=iobs, latmax=latmax, latmin=5, spin=spin, emin=emin_model)
		call_model(fitsfile=fitsfile, iobs=iobs, latmax=latmax, latmin=5,
		spin=spin, emin=emin_model, emax=emax_model)

		filename = 'polarimetry.dat'

		@@ -357,9 +412,10 @@ def write_table_value(latmax=None, ktbb=None, kte=None, tau=None, u_obs=None, sp
		#Here read the I,Q,U spectra of the C program
		#=====================================================================================

		I_flux, Q_flux, U_flux = read_spectra(array_ene=ene_center, filename=filename)
		F_flux, Q_flux, U_flux = read_spectra(array_ene=ene_center, filename=filename)

		#=====================================================================================

		if xspec is True:

		for tt in range(3):
		@@ -371,7 +427,7 @@ def write_table_value(latmax=None, ktbb=None, kte=None, tau=None, u_obs=None, sp


		if tt == 0:
		testspec.setFlux(I_flux)
		testspec.setFlux(F_flux)
		elif tt == 1:
		testspec.setFlux(Q_flux)
		else:
		@@ -380,9 +436,8 @@ def write_table_value(latmax=None, ktbb=None, kte=None, tau=None, u_obs=None, sp
		TableXspec.pushSpectrum(testspec)


		return C_model_kev, PD_filtered


		return matrix_data
		#return C_model_kev, PD_filtered, C_matrix, Q_matrix

		if module_heasp==True:
		TableXspec.setEnergies(logbins_model)
		@@ -418,55 +473,76 @@ if runtable=='n':
		print('with command mcfile=filename')
		exit()

		if ktbb is None:
		print('\nError: ktbb must be defined in order to compute the normalized specific intensity of the 2D fits file')
		print('Set the temperature with ktbb=value')
		os.sys.exit(1)

		if param is None:
		print('Please select the parameter to be plotted with param=name')
		print('Allowed values are C, Q or pd')
		exit(1)



		allowed_params=['C', 'pd']

		if param not in allowed_params:
		print('Unknow param value')
		print('Allowed values are ', allowed_params)
		exit(1)


		#=========================================================================
		#Get parameters of the simulation from filename
		# =========================================================================

		ktbb,kte,tau=filename2params(filename=mcfile)

		if ktbb is None:
		print('Error: ktbb can not be determined from the name of file %s' % (mcfile))
		print('Actually it must be defined in order to compute the normalized '
		'specific intensity of the 2D fits file')
		print('Set the temperature with ktbb=value')
		exit()


		pdf_name = os.path.basename(mcfile.replace(".dat", ".pdf"))
		pdf_name = pdf_name.replace("test_BL", param)

		C_model_kev, PD_filtered = write_table_value(latmax=None, ktbb=ktbb, kte=kte,
		tau=tau, u_obs=None, mc_file=mcfile,
		xspec=False, fitsfile=fitsfile, emin_model=emin_model,
		rebin=rebin, factor=factor)
		#====================================================================================
		#Call the function returning the best fit model and the arrays of data
		#====================================================================================

		result = write_table_value(latmax=None, ktbb=ktbb, kte=kte,tau=tau, u_obs=None, mc_file=mcfile, xspec=False,
		fitsfile=fitsfile, rebin=rebin, factor=factor, energy_center=ene_center,
		delta_energy=delta_energy)

		#print(result.__dict__.keys())

		C_matrix, Q_matrix = read_mc_simulations(mc_file=mcfile, n_i=20, n_j=256)

		plt.figure(figsize=(10, 5))
		#C_model_kev=result[0]
		#PD_filtered=result[1]
		#C_matrix=result[2]
		#Q_matrix=result[3]

		#==========================================================================

		if viewobs is None:
		u_idx=random.randint(0, len(u_centers)-1)
		viewobs=np.arccos(u_centers[u_idx])/np.pi*180

		else:
		u_idx = np.searchsorted(u_bin_edges, cos_viewobs, side='right') - 1
		u_val=u_centers[u_idx]

		# ==========================================================================
		mc_string='__IQ_ADMP'

		if mc_string in pdf_name:
		pdf_name = pdf_name.replace(mc_string, "_iobs"+str(int(viewobs))+mc_string)
		else:
		pdf_name = pdf_name.replace(".pdf", "_iobs" + str(int(viewobs)) + ".pdf")

		#print(u_bin_edges)
		#print(cos_viewobs)


		if pdf is True:
		print('PDF File')
		print(pdf_name)

		@@ -479,61 +555,42 @@ if runtable=='n':
		print('Angle %4.3f deg' % (viewobs))
		# print('==========================')

		#=============================================================================================

		if plot=='y':

		for ii in range(u_idx, u_idx+1):

		if param == 'Q':

		plot_data(x_value=ene_center, y_mc=Q_matrix[ii], y_smooth=Q_filtered[ii],
		title=input_file, ymin=-5500, ymax=5000, logscale=True, pdf_name=pdf_name, pdf=pdf)
		title=input_file, ymin=-5500, ymax=5000, logscale=True,
		pdf_name=pdf_name, pdf=pdf)

		elif param == 'pd':

		safe_Cmatrix = np.where(C_matrix == 0, 1e-20, C_matrix)


		#===================================================================
		#Case of rebin
		#===================================================================

		if rebin==True:

		safe_Cmatrix, ene_center_new, delta_energy_new = fixed_rebin(energy=ene_center, delta_energy=delta_energy,
		Y_array=safe_Cmatrix, factor=factor)
		Q_matrix, ene_center_new, delta_energy_new = fixed_rebin(energy=ene_center, delta_energy=delta_energy,
		Y_array=Q_matrix, factor=factor)

		PD_function = np.interp(ene_center_new, ene_center, PD_filtered[ii])

		ene_center = ene_center_new
		delta_energy = delta_energy_new

		else:

		PD_function= PD_filtered[ii]

		# ===================================================================

		PD_matrix = -Q_matrix / safe_Cmatrix
		error = np.sqrt(2 / safe_Cmatrix[ii])

		plot_data(x_value=ene_center, y_mc=PD_matrix[ii] * 100, y_smooth=PD_function * 100, ymin=ymin,
		ymax=ymax, emin=0.3, emax=emax, y_err=error * 100, y_label='Polarization degree (%)', pdf_name=pdf_name,
		delta_energy=delta_energy, param=param, pdf=pdf)

		Q_interp = np.interp(result.energy_center_rebin, result.energy_center, result.Q_model[ii] *100)

		plot_data(energy_center=result.energy_center_rebin, y_mc=result.PD_matrix_rebin[ii] * 100,
		y_smooth=Q_interp,
		ymin=ymin, ymax=ymax, emin=0.3, emax=30, y_err=result.PD_matrix_rebin_err[ii] * 100,
		y_label='Polarization degree (%)', pdf_name=pdf_name,
		delta_energy=result.delta_energy_rebin, param=param, pdf=pdf)

		elif param == 'C':

		#The division of counts by delta_energy is performed in the plotting routine

		plot_data(x_value=ene_center, y_mc=C_matrix[ii], y_smooth=C_model_kev[ii],
		logscale=True, ymin=0.1, ymax=1e8, y_err=np.sqrt(C_matrix[ii]),
		plot_data(energy_center=ene_center, y_mc=result.C_matrix[ii],
		y_smooth=result.C_model_kev[ii],ymin=0.1, ymax=1e8, y_err=np.sqrt(result.C_matrix[ii]),
		delta_energy=delta_energy, pdf_name=pdf_name,
		param=param, emax=emax)
		param=param, emax=100)


		#Script which creates the FITS file of the simulation
		simul2fits(ene_center=ene_center, delta_ene=delta_energy, counts=C_matrix[ii], err_counts=np.sqrt(C_matrix[ii]), mc_file=mcfile)
		simul2fits(ene_center=ene_center, delta_ene=delta_energy, counts=result.C_matrix[ii],
		err_counts=np.sqrt(result.C_matrix[ii]), mc_file=mcfile)


		else:
		@@ -556,9 +613,10 @@ if runtable=='n':

		print('=================================================================================')
		print('\nIf you want to calculate the I,Q,U spectra for the BL around NS now from the prompt you can type:\n')
		print('bl_pol geom=bl enedep=y enefile=%s iobs=value latmax=value spin=value [emin=emin]' % (fitsfile))
		print('\nDefault minimum energy from the 2D fits file is %3.1f keV' % (emin))
		print('If you want to set lower energy boundary for the back-tracing algorithm set the emin parameter')
		print('bl_pol geom=bl enedep=y enefile=%s iobs=value latmax=value spin=value [emin=emin] [emax=emax]' % (fitsfile))
		print('\nDefault energy interval of the %s file is %3.1f - %3.1f keV' % (fitsfile, emin, emax))

		#print('If you want to set lower energy boundary for the back-tracing algorithm set the emin parameter')


		#====================================================================
		@@ -588,7 +646,7 @@ table_init(init_val=lat_init, array_values=latmax_array, name='latmax', table=Ta
		#====================================================================

		ktbb_array = np.linspace(0.5,3,6)
		ktbb_array = np.arange(1, 3.5, 0.5) #Last value is excluded by arange
		ktbb_array = np.arange(0.5, 3.5, 0.5) #Last value is excluded by arange

		formatted_ktbb_val=[f"{val:.1f}" for val in ktbb_array]
		ktbb_init=2
		@@ -596,7 +654,7 @@ ktbb_init=2
		table_init(init_val=ktbb_init, array_values=ktbb_array, name='kTbb', table=TableXspec)

		#kte_array = np.linspace(1.5,5,8)
		kte_array=np.array([2,2.5,3,3.5,4,4.5,5,6,7])
		kte_array=np.array([1,2,3,4,5,6,7])
		formatted_kte_val=[f"{val:.1f}" for val in kte_array]
		kte_init=3

		@@ -608,7 +666,7 @@ table_init(init_val=kte_init, array_values=kte_array, name='kTe', table=TableXsp

		tau_array = np.linspace(0.5,4, 8)

		tau_array=np.array([1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 6])
		tau_array=np.array([0.5,1,2,3,4,5,6,7])

		formatted_tau_val=[f"{val:.1f}" for val in tau_array]
		tau_init=4
		@@ -616,7 +674,6 @@ tau_init=4
		table_init(init_val=tau_init, array_values=tau_array, name='tau', table=TableXspec)



		#====================================================================
		#cos(i_obs)
		#====================================================================
		@@ -666,14 +723,14 @@ for latmax in latmax_array:

		for flag in flag_array:

		filename='test_BL_slab_Z0seed_grid__'+ktbb+'_'+kte+'_'+tau+'__IQ_ADMP.dat'
		#filename='test_BL_slab_Z0seed_grid__'+ktbb+'_'+kte+'_'+tau+'__IQ_ADMP.dat'
		filename='test_iasfBO_'+ktbb+'-'+kte+'-'+tau+'_slab_seedZ0_pol0_IQ_dump.dat'

		if flag==0:
		file_path = '../unpolarized_seed/'+filename
		else:
		file_path = '../polarized_seed/'+ filename


		paramvalues = np.array([float(latmax), float(ktbb), float(kte), float(tau), u_obs, flag])

		#paramvalues = np.array([float(latmax), float(ktbb), float(kte), float(tau), u_obs, spin, flag])
		@@ -685,8 +742,10 @@ for latmax in latmax_array:
		print('Check the setup of the variable FreeParams in the script')
		os.sys.exit(1)

		write_table_value(latmax=latmax, ktbb=float(ktbb), kte=float(kte), tau=float(tau), u_obs=u_obs, spin=spin, mc_file=file_path,
		xspec=True, paramvalues=paramvalues, emin_model=emin_model, rebin=rebin)
		write_table_value(latmax=latmax, ktbb=float(ktbb), kte=float(kte), tau=float(tau),
		u_obs=u_obs, spin=spin, mc_file=file_path,
		xspec=True, paramvalues=paramvalues,
		emin_model=emin_model, emax_model=emax_model, rebin=rebin, factor=factor)


		# now write out the table.

build_model/model_utilities.py

+688 −129

File changed.

Preview size limit exceeded, changes collapsed.

build_model/rebin_spectrum.py

+69 −1

Original line number	Diff line number	Diff line
		import numpy as np


		def log_rebin(Y_array=None, energy=None, delta_energy=None, N_bins=None):
		"""
		Ribinaggio logaritmico di Y_array usando energy e delta_energy.
		Supporta Y_array 1D o 2D (più spettri).

		Parametri:
		Y_array : np.ndarray
		Array dei dati da ribinnare (1D o 2D)
		energy : np.ndarray
		Centri dei bin originali
		delta_energy : np.ndarray
		Larghezze dei bin originali
		N_bins : int
		Numero di bin logaritmici da creare

		Ritorna:
		Y_new : np.ndarray
		Array ribinnato
		energy_new : np.ndarray
		Centri dei nuovi bin logaritmici
		delta_energy_new : np.ndarray
		Larghezze dei nuovi bin
		"""

		if Y_array.ndim == 1:
		Y_array = Y_array[:, np.newaxis]

		Nview_obs, Nbins_ene = Y_array.shape

		# Calcolo dei bordi originali dai centri e delta_energy
		edges_min = energy - delta_energy / 2
		edges_max = energy + delta_energy / 2

		E_min = edges_min[0]
		E_max = edges_max[-1]

		# Bordi dei nuovi bin logaritmici
		edges_new = np.logspace(np.log10(E_min), np.log10(E_max), N_bins + 1)

		# Array output
		Y_new = np.zeros((Nview_obs, N_bins))
		energy_new = np.zeros(N_bins)
		delta_energy_new = np.zeros(N_bins)

		for ii in range(Nview_obs):
		for jj in range(N_bins):
		E_min_bin = edges_new[jj]
		E_max_bin = edges_new[jj + 1]

		# Trova i bin originali che rientrano nel nuovo bin
		idx = (edges_min >= E_min_bin) & (edges_max < E_max_bin)

		# Somma dei valori
		Y_new[ii, jj] = Y_array[ii, idx].sum()

		# Centro e larghezza dei nuovi bin
		if np.any(idx):
		energy_new[jj] = 0.5 * (edges_min[idx][0] + edges_max[idx][-1])
		else:
		energy_new[jj] = 0.5 * (E_min_bin + E_max_bin)
		delta_energy_new[jj] = E_max_bin - E_min_bin

		if Nview_obs == 1:
		Y_new = Y_new.flatten()

		return Y_new, energy_new, delta_energy_new

		def fixed_rebin(energy=None, delta_energy=None, Y_array=None, factor=None):
		"""
		Ribinaggio fisso dei dati in energia.
		@@ -14,6 +82,7 @@ def fixed_rebin(energy=None, delta_energy=None, Y_array=None, factor=None):

		Nview_obs, Nbins_ene = Y_array.shape # N = numero di bin, M = numero di spettri


		Nbins_ene_new=int(Nbins_ene/factor)

		# Bordi originali
		@@ -63,7 +132,6 @@ def fixed_rebin(energy=None, delta_energy=None, Y_array=None, factor=None):

		return Y_new, energy_new, delta_energy_new


		def adaptive_rebin_sn(energy, Y, Yerr, target_sn):
		"""
		Rebin adattivo: unisce canali finché il S/N del bin ≥ target_sn.