Defined a new parameter at the shell prompt for grouping the Q spectra (f4762d9c) · Commits · polarization / boundary_layer

build_model/REAME.txt

0 → 100644

+17 −0

Original line number	Diff line number	Diff line
		With current version of the code the left boundary of the energy channels
		of the MC simulation is 0.3 keV

		1) To create the FITS table model with lower boundary, use the following command:

		make_model.py runtable=y tabname=blcomp_new.fits emin_model=0.1

		To rebin the Q spectra for both fitting and plotting purposes, add to the command line
		rebin=y

		2) To show single plot with fit

		make_model.py runtable=n mcfile=../unpolarized_seed/test_BL_slab_Z0seed_grid__2.5_7.0_3.0__IQ_ADMP.dat ktbb=2 param=pd ymin=-8 ymax=5 emax=60 viewobs=70 emax=30 [pdf=y]

		the latest parameter allows to create the PDF figure from the plot

build_model/comptb.py

100644 → 100755

+0 −0

File mode changed from 100644 to 100755.

View file

build_model/comptb_numba.py

100644 → 100755

+0 −0

File mode changed from 100644 to 100755.

View file

build_model/make_model.py

+76 −16

Original line number	Diff line number	Diff line
		@@ -6,6 +6,7 @@ import random
		from scipy.optimize import curve_fit
		from comptb import *
		from model_utilities import *
		from rebin_spectrum import *


		#=======================================================================
		@@ -66,6 +67,15 @@ for inputval in sys.argv:
		if inputval.startswith('emin_model') == True:
		emin_model = float(inputval[inputval.rfind('=') + 1:])

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

		if inputval.startswith('factor') == True:
		factor = int(inputval[inputval.rfind('=') + 1:])

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

		try:
		runtable
		except:
		@@ -123,6 +133,35 @@ try:
		except:
		timeon='n'

		try:
		rebin

		if rebin=='y':
		rebin=True
		elif rebin=='n':
		rebin=False
		else:
		print('rebin parameter not set correctly')
		exit()
		except:
		rebin=False

		try:
		factor
		except:
		factor=4

		try:
		pdf

		if pdf=='y':
		pdf=True
		else:
		pdf=False

		except:
		pdf=False


		#=====================================================================================
		#Energy bins
		@@ -245,7 +284,7 @@ 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):
		xspec=None, fitsfile=None, paramvalues=None, emin_model=None, rebin=None, factor=None):

		print('==============================================')
		print('Reading file: %s' % (mc_file))
		@@ -258,7 +297,6 @@ def write_table_value(latmax=None, ktbb=None, kte=None, tau=None, u_obs=None, sp
		C_model_kev = smooth_matrix_spectrum(specpol_array=C_matrix, energy_center=ene_center,
		delta_energy=delta_energy, ktbb=ktbb, kte=kte, tau=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

		@@ -269,7 +307,8 @@ def write_table_value(latmax=None, ktbb=None, kte=None, tau=None, u_obs=None, sp
		# ===============================================================================

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


		#print(C_matrix.shape, Q_matrix.shape, C_model_kev.shape, PD_filtered.shape)
		@@ -407,12 +446,11 @@ if runtable=='n':

		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)
		xspec=False, fitsfile=fitsfile, emin_model=emin_model,
		rebin=rebin, factor=factor)

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



		plt.figure(figsize=(10, 5))

		if viewobs is None:
		@@ -437,34 +475,56 @@ if runtable=='n':
		# print(" %d %4.3f %4.3f" % (ii, u_centers[ii], np.arccos(u_centers[ii])/np.pi*180))

		# print('==========================')
		print("Selected window ", )
		print("Selected window %d" % (u_idx))
		print('Angle %4.3f deg' % (viewobs))
		# print('==========================')

		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)
		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)

		PD_matrix = -Q_matrix / safe_Cmatrix

		#===================================================================
		#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_filtered[ii] * 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)
		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)

		elif param == 'C':

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

		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]),
		@@ -626,7 +686,7 @@ for latmax in latmax_array:
		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)
		xspec=True, paramvalues=paramvalues, emin_model=emin_model, rebin=rebin)


		# now write out the table.

build_model/model_utilities.py

+90 −29

Original line number	Diff line number	Diff line
		@@ -2,12 +2,10 @@ import math
		from fileinput import filename
		from sndhdr import test_au
		import numpy as np
		import inspect
		import os
		import sys
		import re
		import warnings
		import random
		from astropy.io import fits
		from scipy.signal import savgol_filter
		from scipy.optimize import curve_fit,minimize_scalar
		@@ -23,6 +21,9 @@ from lmfit import Model
		from matplotlib.ticker import AutoMinorLocator
		from comptb_numba import *
		from fit_functions import *
		from rebin_spectrum import *
		from matplotlib.ticker import MultipleLocator


		warnings.filterwarnings("ignore", category=RuntimeWarning)

		@@ -505,7 +506,11 @@ def fit_log_data(energy_center=None, array_y=None, y_err=None, deg=None, emin=N

		mask = (energy_center > emin) & (energy_center < emax) & (array_y > 0)



		if y_err is None:


		y_err = np.sqrt(np.abs(array_y))
		safe_y_err = np.where(y_err == 0, 1e8, y_err)

		@@ -630,7 +635,7 @@ def angular_integration(u_array=None, Y_array=None):

		#======================================================================
		def plot_data(x_value=None, y_mc=None, y_smooth=None, title=None, ymin=None, ymax=None, emin=None, emax=None,
		logscale=None, y_err=None, y_label=None, delta_energy=None, pdf_name=None, param=None):
		logscale=None, y_err=None, y_label=None, delta_energy=None, pdf_name=None, param=None, pdf=None):


		mask= (x_value > 0.3) & (x_value < emax)
		@@ -651,8 +656,6 @@ def plot_data(x_value=None, y_mc=None, y_smooth=None, title=None, ymin=None, yma
		else:
		ymin=1.1*ymin_mc



		if ymax is None:
		if ymax_mc > 0:
		ymax = 1.2 * ymax_mc
		@@ -705,6 +708,7 @@ def plot_data(x_value=None, y_mc=None, y_smooth=None, title=None, ymin=None, yma
		ax1.set_xscale('log')
		ax1.set_ylim(bottom=ymin, top=ymax)


		if param=='C':
		ax1.set_yscale('log')
		ax1.set_ylabel(r"Counts keV$^{-1}$", fontsize=14)
		@@ -713,9 +717,12 @@ def plot_data(x_value=None, y_mc=None, y_smooth=None, title=None, ymin=None, yma


		elif param=='pd':

		ax1.set_ylabel(r"Q (%)", fontsize=14)
		ax1.yaxis.set_major_formatter(ticker.FormatStrFormatter('%d'))



		ax1.set_xlim(left=0.32, right=emax)
		ax1.set_title(title)

		@@ -725,7 +732,6 @@ def plot_data(x_value=None, y_mc=None, y_smooth=None, title=None, ymin=None, yma
		label.set_fontsize(11)



		ax2.set_xscale('log')
		ax2.axhline(0, color='gray', lw=1)
		ax2.set_ylim(bottom=-4, top=4)
		@@ -769,19 +775,28 @@ def plot_data(x_value=None, y_mc=None, y_smooth=None, title=None, ymin=None, yma
		data=y_mc
		data_err=y_err



		ax1.plot(x_value, model, marker='o', markersize=1, label='', linewidth=2)
		ax1.errorbar(x_value, data, yerr=data_err, xerr=delta_energy/2, fmt='s', markersize=2,
		linewidth=1.5, capsize=0, label='Monte Carlo', color='orange')
		linewidth=1.5, capsize=0, label='Monte Carlo', color='black')

		resid= (model-data)/data_err
		resid_err=np.ones(len(x_value))

		ax2.errorbar(x_value, resid, yerr=resid_err, xerr=delta_energy/2, fmt='s', markersize=2,
		linewidth=1.5, capsize=0, label='Monte Carlo', color='blue')
		linewidth=1.5, capsize=0, label='Monte Carlo', color='black')


		if param=='pd':
		y_ticks = np.arange(int(np.floor(ymin)), int(np.ceil(ymax)) + 1, 1)
		y_ticks = np.arange(int(np.floor(ymin)), int(np.ceil(ymax)) + 2, 2)
		#y_ticks = np.arange(-10, 10 + 1, 2)

		minor_ticks = np.arange(int(np.floor(ymin)), int(np.ceil(ymax)) , 1)
		ax1.set_yticks(minor_ticks, minor=True)
		ax1.minorticks_on()



		else:
		exp_min = int(np.floor(np.log10(ymin)))
		@@ -803,11 +818,11 @@ def plot_data(x_value=None, y_mc=None, y_smooth=None, title=None, ymin=None, yma
		plt.errorbar(x_value, y_mc / delta_energy, yerr=y_err / delta_energy, fmt='s', markersize=2,
		linewidth=1.5, capsize=2, label='Monte Carlo', color='orange')

		if pdf_name is not None:
		if pdf_name is not None and pdf is True:
		plt.savefig(pdf_name)


		#plt.tight_layout()
		plt.tight_layout()
		#plt.legend()

		#==========================================================================
		@@ -846,7 +861,7 @@ def read_spectra(array_ene=None, filename=None):
		return I_flux_photons,Q_flux_photons, U_flux_photons

		#====================================================================
		def read_mc_simulations(mc_file=None, n_i=20, n_j=256):
		def read_mc_simulations(mc_file=None, n_i=20, n_j=256, energy_center=None, delta_energy=None, rebin=None):

		if mc_file is not None:
		if not os.path.exists(mc_file):
		@@ -893,6 +908,14 @@ def read_mc_simulations(mc_file=None, n_i=20, n_j=256):
		Q_grouped = Q.reshape(10, 2, 256).sum(axis=1)
		C_grouped = C.reshape(10, 2, 256).sum(axis=1)

		if rebin is True:
		if energy_center is not None and delta_energy is not None:
		print('CIAO')

		else:
		print('Error: either the array of energy centers or the array of delta_energies is not define')
		exit(1)


		print('Returning grouped C and Q spectra')
		return C_grouped, Q_grouped
		@@ -1024,11 +1047,9 @@ def smooth_matrix_spectrum(specpol_array=None, energy_center=None, delta_energy=
		#Note that the array with fit is divided by delta energy, so its units are counts/keV
		#========================================================================================


		smooth_specpol_array[ii] = fit_poly_cpl(energy_center=energy_center, delta_energy=delta_energy,
		array_y=specpol_array[ii], ktbb=ktbb, kte=kte, tau=tau)


		if np.all(smooth_specpol_array[ii] == False):
		print('Failed to converge')

		@@ -1046,21 +1067,60 @@ def residuals_const_piecewise_parabola(params, x, y, yerr):
		return (model - y) / yerr


		def smooth_matrix_pd(energy_center=None, Q_matrix=None, C_matrix=None, emin=None, emax=20):
		def smooth_matrix_pd(energy_center=None, delta_energy=None, Q_matrix=None, C_matrix=None, emin=None,
		emax=20, rebin=None, factor=None):

		row, col = Q_matrix.shape
		smooth_pd_matrix = np.zeros_like(Q_matrix)
		if rebin is True and factor is None:
		print('Error in routine smooth_matrix_pd: rebinning energy factor is not define')
		exit(1)

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

		#=======================================================
		#Rebin the C and Q matrixes along the energy axis
		#=======================================================

		if rebin==True:

		Q_matrix_rebin, energy_new, delta_energy_new = fixed_rebin(energy=energy_center, delta_energy=delta_energy,
		Y_array=Q_matrix, factor=factor)

		safe_Cmatrix_rebin, energy_center_new, delta_energy_new = fixed_rebin(energy=energy_center, delta_energy=delta_energy,
		Y_array=safe_Cmatrix, factor=factor)

		energy_center_local = energy_center_new
		delta_energy_local = delta_energy_new

		pd_matrix = Q_matrix_rebin / safe_Cmatrix_rebin
		pd_error = np.sqrt(2 / safe_Cmatrix_rebin)

		print("Fitting the rebinned Q spectrum")


		else:

		energy_center_local=energy_center
		delta_energy_local = delta_energy

		pd_matrix = Q_matrix / safe_Cmatrix
		pd_error = np.sqrt(2 / safe_Cmatrix)

		#================================================================
		#The final smoothed PD matrix for each viewing angle will have
		#the same number of energy channels as the original
		#================================================================

		row, col = Q_matrix.shape
		smooth_pd_matrix = np.zeros_like(Q_matrix)

		#print("Dimensioni Q_array rebinnato:", Q_matrix.shape)
		#print("Dimensioni C_array rebinnato:", safe_Cmatrix.shape)

		#For the error on Q, see Eq. (23) of Kislat et al. 2015
		#https://arxiv.org/pdf/1409.6214

		error = np.sqrt(2 / safe_Cmatrix)

		mask = (energy_center > emin) & (energy_center <= emax)

		mask = (energy_center_local-delta_energy_local > emin) & (energy_center_local+delta_energy_local <= emax)

		for ii in range(row):

		@@ -1069,9 +1129,9 @@ def smooth_matrix_pd(energy_center=None, Q_matrix=None, C_matrix=None, emin=None

		for deg in degrees:

		x_value = energy_center[mask]
		x_value = energy_center_local[mask]
		y_value = pd_matrix[ii][mask]
		y_err = error[ii][mask]
		y_err = pd_error[ii][mask]

		n = len(x_value)

		@@ -1084,8 +1144,6 @@ def smooth_matrix_pd(energy_center=None, Q_matrix=None, C_matrix=None, emin=None

		results.append({'deg': deg, 'chi2': chi2, 'BIC': BIC, 'poly': poly})



		#=================================================================
		#Stores results for increaing BIC
		#=================================================================
		@@ -1095,17 +1153,20 @@ def smooth_matrix_pd(energy_center=None, Q_matrix=None, C_matrix=None, emin=None
		#for r in results_sorted:
		# print(f"deg: {r['deg']}, chi2: {r['chi2']}, BIC: {r['BIC']}")


		# Modello migliore
		best_fit = results_sorted[0]
		chi2_poly = best_fit['chi2']


		smooth_pd_matrix[ii] = best_fit['poly'](energy_center)
		if rebin==True:

		smooth_pd_matrix[ii] = np.interp(energy_center, energy_center_local, best_fit['poly'](energy_center_local))

		else:
		smooth_pd_matrix[ii] = best_fit['poly'](energy_center_local)

		return smooth_pd_matrix

		return smooth_pd_matrix

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