Written new routine for adaptive grouping of the polarization spectra (729ccc12) · Commits · polarization / boundary_layer

build_model/REAME.txt

+9 −4

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:
		By default the Q spectra are rebinned by a factor equal to 8 with respect to the original 256 channels,
		to improve the S/N.

		To see the original unbinned Q spectra at the prompt type rebin=n, while to change the binning
		type factor=M where the latter must be a power of 2.


		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]
		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 viewobs=70 emax=30 [pdf=y]

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

build_model/make_model.py

+81 −42

Original line number	Diff line number	Diff line
		@@ -25,61 +25,67 @@ except ImportError:
		#====================================================================
		for inputval in sys.argv:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

		try:
		@@ -177,7 +183,7 @@ except:
		#which will be used by the C program bl_pol to perform extrapolation
		#====================================================================================

		energy_bins=257

		try:
		emin
		except:
		@@ -198,18 +204,29 @@ try:
		except:
		emax_model=emax

		try:
		emin_plot
		except:
		emin_plot=emin

		try:
		emax_plot
		except:
		emax_plot=emax

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

		energy_bins=259

		logbins = np.logspace(np.log10(emin), np.log10(emax), num=energy_bins)

		ene_center=np.asarray([round((logbins[ii]+logbins[ii+1])/2,4) for ii in range(energy_bins-1)])
		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_model), num=energy_bins-2)

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

		ene_center = ene_center[1:-1]
		delta_energy = delta_energy[1:-1]

		#============================================================================
		#Set table descriptors and the energy array
		@@ -278,8 +295,8 @@ def simul2fits(ene_center=None, delta_ene=None, counts=None, err_counts=None, mc
		#========================================================================
		def table_init(init_val=None, array_values=None, name=None, table=None):

		min_val=array_values[0]
		max_val=array_values[-1]
		min_val=float(array_values[0])
		max_val=float(array_values[-1])

		testpar = tableParameter()
		testpar.setName(name)
		@@ -288,6 +305,8 @@ def table_init(init_val=None, array_values=None, name=None, table=None):
		testpar.setInitialValue(init_val)
		testpar.setDelta(0.1)



		testpar.setMinimum(min_val)
		testpar.setBottom(min_val)

		@@ -318,10 +337,11 @@ def write_table_value(latmax=None, ktbb=None, kte=None, tau=None, u_obs=None, sp
		print('Parameters value', ktbb, kte, tau)


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


		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)

		@@ -355,6 +375,8 @@ def write_table_value(latmax=None, ktbb=None, kte=None, tau=None, u_obs=None, sp
		#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,
		@@ -364,6 +386,7 @@ def write_table_value(latmax=None, ktbb=None, kte=None, tau=None, u_obs=None, sp

		matrix_data.Q_model = PD_filtered


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

		@@ -379,6 +402,7 @@ def write_table_value(latmax=None, ktbb=None, kte=None, tau=None, u_obs=None, sp
		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'

		@@ -388,6 +412,8 @@ def write_table_value(latmax=None, ktbb=None, kte=None, tau=None, u_obs=None, sp
		print('Error while creating the fits file %s ' % (fitsfile))
		os.sys.exit(1)



		if spin is None:
		spin=500

		@@ -416,6 +442,7 @@ def write_table_value(latmax=None, ktbb=None, kte=None, tau=None, u_obs=None, sp

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


		if xspec is True:

		for tt in range(3):
		@@ -512,7 +539,7 @@ if runtable=='n':
		fitsfile=fitsfile, rebin=rebin, factor=factor, energy_center=ene_center,
		delta_energy=delta_energy)

		#print(result.__dict__.keys())



		#C_model_kev=result[0]
		@@ -531,12 +558,12 @@ if runtable=='n':
		u_val=u_centers[u_idx]

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

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

		#print(u_bin_edges)
		#print(cos_viewobs)
		@@ -570,22 +597,27 @@ if runtable=='n':
		elif param == 'pd':


		Q_interp = np.interp(result.energy_center_rebin, result.energy_center, result.Q_model[ii] *100)
		#Q_interp = np.interp(result.energy_center_rebin[ii], 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,
		Q_model = result.Q_model[ii] * 100

		plot_data(energy_center_data=result.energy_center_rebin[ii], y_mc=result.PD_matrix_rebin[ii] * 100,
		y_smooth=Q_model, energy_center_model=ene_center,
		ymin=ymin, ymax=ymax, emin=0.3, emax=35, 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)
		delta_energy_data=result.delta_energy_rebin[ii], param=param,
		pdf=pdf)

		elif param == 'C':


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

		plot_data(energy_center=ene_center, y_mc=result.C_matrix[ii],

		plot_data(energy_center_data=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=100)
		delta_energy_data=delta_energy, pdf_name=pdf_name,
		param=param, pdf=pdf, emax=emax_plot)


		#Script which creates the FITS file of the simulation
		@@ -633,6 +665,8 @@ if runtable=='n':
		#Here define the grid values for kTbb, kTe, tau and i_obs
		#====================================================================



		#Latmax
		latmax_array = np.array([10., 30., 50., 70., 90.])

		@@ -646,7 +680,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(0.5, 3.5, 0.5) #Last value is excluded by arange
		ktbb_array = np.arange(1.0, 3.5, 0.5) #Last value is excluded by arange

		formatted_ktbb_val=[f"{val:.1f}" for val in ktbb_array]
		ktbb_init=2
		@@ -666,14 +700,13 @@ 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([0.5,1,2,3,4,5,6,7])
		tau_array=np.array([1,2,3,4,5,6,7])

		formatted_tau_val=[f"{val:.1f}" for val in tau_array]
		tau_init=4

		table_init(init_val=tau_init, array_values=tau_array, name='tau', table=TableXspec)


		#====================================================================
		#cos(i_obs)
		#====================================================================
		@@ -691,7 +724,7 @@ spin_array=np.linspace(0,700, 5)
		spin_init=300
		spin=500

		spin_array=np.array([1000])
		spin_array=np.array([500])
		#table_init(init_val=spin_init, array_values=spin_array, name='spin', table=TableXspec)

		#====================================================================
		@@ -724,13 +757,17 @@ for latmax in latmax_array:
		for flag in flag_array:

		#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:
		filename = 'test_iasfBO_' + ktbb + '-' + kte + '-' + tau + '_slab_seedZ0_pol0_IQ_dump.dat'
		file_path = '../unpolarized_seed/'+filename
		else:

		filename = 'test_iasfBO_' + ktbb + '-' + kte + '-' + tau + '_slab_seedZ0_pol42_IQ_dump.dat'
		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])
		@@ -745,7 +782,9 @@ for latmax in latmax_array:
		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)
		emin_model=emin_model,
		emax_model=emax_model,
		rebin=rebin, factor=factor, energy_center=ene_center,delta_energy=delta_energy)


		# now write out the table.

build_model/model_utilities.py

+420 −82

File changed.

Preview size limit exceeded, changes collapsed.

build_model/rebin_spectrum.py

+101 −25

Original line number	Diff line number	Diff line
		@@ -132,28 +132,104 @@ 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.
		Restituisce: energy_new, Y_new, Yerr_new, snr_new, bin_edges
		"""
		N = len(Y)
		i = 0
		energy_new, Y_new, Yerr_new, snr_new, bin_edges = [], [], [], [], []
		while i < N:
		sumY, sumYerr2 = 0.0, 0.0
		j = i
		while j < N:
		sumY += Y[j]
		sumYerr2 += Yerr[j]**2
		cur_snr = sumY / np.sqrt(sumYerr2) if sumYerr2 > 0 else 0.0
		j += 1
		if cur_snr >= target_sn:

		#================================================================================================
		#Adaptive rebinning: for each Q energy channel require to have target_sn value
		#assuming P_deg=1/% and error on P_deg given by Sqrt[2/N]
		#P/Sqrt[2/N]=SNR so that N=2SNR2/P_deg*2 where N is the number of counts in the given bin
		#if the minimum number of counts is not achieved, group a maximum of 8 bins per new channel
		#================================================================================================

		def adaptive_rebin_sn(C_matrix=None, Q_matrix=None, energy=None, delta_energy=None, target_sn=3):


		P_deg=0.01
		Min_counts=2 * target_sn2/P_deg2

		if energy is None:
		print('Undefined energy array in routine adaptive_rebin_sn')
		exit(1)


		Nview_obs, Nbins_ene = C_matrix.shape

		Q_rebinned_all = []
		C_rebinned_all = []
		energy_rebinned_all = []
		delta_energy_rebinned_all = []

		if np.any(C_matrix == 0):
		print("L'array C_matrix contiene almeno uno zero")
		else:
		print("Nessuno zero in C_matrix")


		for ii in range(Nview_obs):

		q_row = Q_matrix[ii, :]
		c_row = C_matrix[ii, :]

		Q_rebinned_row = []
		C_rebinned_row = []
		energy_rebinned_row = []
		delta_energy_rebinned_row = []

		start_bin = 0

		while start_bin < Nbins_ene:

		cum_c = 0
		bins_to_group = 0

		for jj in range(start_bin, min(start_bin + 8, Nbins_ene)):

		cum_c += c_row[jj]
		bins_to_group += 1

		#print('Angle %d bin %d Counts %d Bins to group %s' % (ii, jj, c_row[jj], bins_to_group))

		if cum_c >= Min_counts or bins_to_group >= 8:

		#print('Achieved', target_sn, bins_to_group)
		break
		energy_new.append(energy[i:j].mean())
		Y_new.append(sumY)
		Yerr_new.append(np.sqrt(sumYerr2))
		snr_new.append(sumY / np.sqrt(sumYerr2) if sumYerr2 > 0 else 0.0)
		bin_edges.append((i, j-1))
		i = j
		return np.array(energy_new), np.array(Y_new), np.array(Yerr_new), np.array(snr_new), bin_edges

		# Accumulo semplice dei valori per C e Q

		sum_Q=np.sum(q_row[start_bin:start_bin + bins_to_group])
		sum_C=np.sum(c_row[start_bin:start_bin + bins_to_group])

		Q_rebinned_row.append(sum_Q)
		C_rebinned_row.append(sum_C)

		# Energia centroide = centro geometrico del nuovo bin
		E_min = energy[start_bin]-delta_energy[start_bin]/2
		E_max = energy[start_bin + bins_to_group - 1]+delta_energy[start_bin + bins_to_group - 1]/2

		E_centroid = (E_min + E_max) / 2
		energy_rebinned_row.append(E_centroid)

		# Larghezza totale del nuovo bin
		delta_energy_rebinned_row.append(np.sum(delta_energy[start_bin:start_bin + bins_to_group]))

		#print('New bin grouped: E_min = {:.3f}, E_max = {:.3f}, E_centroid = {:.3f},'.format(E_min, E_max, E_centroid))

		start_bin += bins_to_group

		Q_rebinned_all.append(np.array(Q_rebinned_row))
		C_rebinned_all.append(np.array(C_rebinned_row))
		energy_rebinned_all.append(np.array(energy_rebinned_row))
		delta_energy_rebinned_all.append(np.array(delta_energy_rebinned_row))

		Q_rebinned_all = np.array(Q_rebinned_all, dtype=object)
		C_rebinned_all = np.array(C_rebinned_all, dtype=object)
		energy_rebinned_all = np.array(energy_rebinned_all, dtype=object)
		delta_energy_rebinned_all = np.array(delta_energy_rebinned_all, dtype=object)



		return Q_rebinned_all, C_rebinned_all, energy_rebinned_all, delta_energy_rebinned_all

main_program.c

+4 −0

Original line number	Diff line number	Diff line
		@@ -635,6 +635,10 @@ int main(int argc, char* argv[])
		/=====================================================================================/
		if (FlagEneDep == TRUE)
		{
		if (kk==0)
		{
		fprintf(dat, "#Emin Emax Flux Q U PA\n");
		}
		fprintf(dat, "%5.4f %5.4f %5.4e %4.3e %4.3e %5.4f\n", array_ene[kk] - array_delta_ene[kk] / 2,
		array_ene[kk] + array_delta_ene[kk] / 2, primary_cap->Ival[ii] * factor,
		primary_cap->Qval[ii] * factor, primary_cap->Uval[ii] * factor, PA_obs / PI * 180);