phangs_compare_seds_isrf.pro
3.77 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
PRO phangs_compare_seds_isrf,isrf_class $
,N_seds=N_seds $
,normalize=normalize $
,weighted=weighted $
,include_herschel_filters=include_herschel_filters
;phangs_compare_seds_isrf,0,N_seds=100L ;test works
;phangs_compare_seds_isrf,0,N_seds=100L,/weighted ;not working because weights are infinite
;phangs_compare_seds_isrf,0,N_seds=100L,/normalize ;not working because weights are infinite
;phangs_compare_seds_isrf,15,N_seds=100L
;phangs_compare_seds_isrf,15,N_seds=1000L,/normalize
define_la_common
dir=!phangs_data_dir+'/ISRF/GRIDS/'
isrf_class_str='_isrfclass'+strtrim(isrf_class,2)
model='DBP90'
table_name=dir+model+'_MuseISRF_JWST_G0_YPAH_YVSG_4Phangs'+isrf_class_str+'.fits'
st=mrdfits(table_name,1,h)
Nst=n_elements(st)
use_filters_names=['F360M','F0770W','F1000W','F1130W','F2100W']
IF keyword_set(include_herschel_filters) THEN BEGIN
use_filters_names=['F360M','F0770W','F1000W','F1130W','F2100W','PACS_BLUE','PACS_GREEN','PACS_RED','PSW','PMW','PLW']
ENDIF
filters=dustem_filter_names2filters(use_filters_names)
Nfilters=n_elements(filters)
ii=100L
Nseds=100L
IF keyword_set(N_seds) THEN Nseds=N_seds
Nparams=3
true_params=fltarr(Nseds,Nparams)
fitted_params=fltarr(Nseds,Nparams)
facts=fltarr(Nseds)
seed=1
FOR j=0L,Nseds-1 DO BEGIN
;seed=1
inn=abs(randomu(seed))
ii=inn*Nst
;print,inn,ii
sed=dustem_initialize_sed(Nfilters)
;sed=fltarr(Nfilters)
tags=tag_names(st)
FOR i=0L,Nfilters-1 DO BEGIN
ind=where(tags EQ 'I'+filters[i],count)
sed[i].filter=filters[i]
sed[i].stokesI=st[ii].(ind[0])
ENDFOR
sed.sigmaii=sed.stokesI*1.e-5
FOR i=0L,Nparams-1 DO BEGIN
true_params[j,i]=st[ii].(i)
ENDFOR
;stop
fit_isrf_class_str='_isrfclass0'
table_name_2=dir+model+'_MuseISRF_JWST_G0_YPAH_YVSG_4Phangs'+fit_isrf_class_str+'.fits'
params=dustem_brute_force_fit(sed,table_name_2,filters,fact=fact,chi2=chi2,normalize=normalize,rchi2=rchi2,show_sed=show_sed $
,params_hit=params_hit,params_uncertainties=params_uncertainties,params_min=params_min,params_max=params_max $
,fixed_parameters_description=fixed_parameters_description,fixed_parameters_values=fixed_parameters_values $
,weighted_params=weighted_params)
print,params,fact
facts[j]=fact
params[1]=params[1]/fact
params[2]=params[2]/fact
weighted_params[1]=weighted_params[1]/fact
weighted_params[2]=weighted_params[2]/fact
;fitted_params[j,*]=params
IF keyword_set(weighted) THEN BEGIN
fitted_params[j,*]=weighted_params
ENDIF ELSE BEGIN
fitted_params[j,*]=params
ENDELSE
ENDFOR
;stop
win=0L
!p.multi=[0,1,3]
xtit='True parameter ('+isrf_class_str+')'
ytit='recovered parameter (Mathis ISRF)'
tit='G0'
tit_str=':'
IF keyword_set(weighted) THEN tit_str=tit_str+' weighted'
IF keyword_set(normalize) THEN tit_str=tit_str+' normalized'
IF keyword_set(include_herschel_filters) THEN tit_str=tit_str+' MUSE+Herschel' ELSE tit_str=tit_str+' MUSE'
;window,win & win=win+1
window,win,xsize=500,ysize=900
xrange=minmax(true_params[*,0])
cgplot,true_params[*,0],fitted_params[*,0],tit='G0'+tit_str,psym=4,/xlog,/ylog,xtit=xtit,ytit=ytit,xrange=xrange,yrange=xrange
cgoplot,10^!x.crange,10^!y.crange,linestyle=2,color='red',thickness=3
;window,win & win=win+1
xrange=minmax(true_params[*,1])
cgplot,true_params[*,1],fitted_params[*,1],tit='Ypah'+tit_str,psym=4,/xlog,/ylog,xtit=xtit,ytit=ytit,xrange=xrange,yrange=xrange
cgoplot,10^!x.crange,10^!y.crange,linestyle=2,color='red',thickness=3
;window,win & win=win+1
xrange=minmax(true_params[*,2])
cgplot,true_params[*,2],fitted_params[*,2],tit='Yvsg'+tit_str,psym=4,/xlog,/ylog,xtit=xtit,ytit=ytit,xrange=xrange,yrange=xrange
cgoplot,10^!x.crange,10^!y.crange,linestyle=2,color='red',thickness=3
stop
cleanplot
window,win & win=win+1
cgplot,facts,true_params[*,1]-fitted_params[*,1],psym=4
END