phangs_compare_seds_isrf.pro
5.47 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
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
PRO phangs_compare_seds_isrf,isrf_class $
,isrf_class_analysis=isrf_class_analysis $
,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 some mismatch
;phangs_compare_seds_isrf,15,N_seds=100L
;phangs_compare_seds_isrf,15,N_seds=1000L,/normalize
;phangs_compare_seds_isrf,15,N_seds=1000L,/normalize,/weighted
;phangs_compare_seds_isrf,15,N_seds=1000L,/normalize,/weighted,/include_herschel_filters
;phangs_compare_seds_isrf,15,N_seds=1000L,/normalize,/weighted,/include_herschel_filters,isrf_class_analysis=14
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)
chi2s=fltarr(Nseds)
seed=1
use_isrf_class_analysis=0L
IF keyword_set(isrf_class_analysis) THEN BEGIN
use_isrf_class_analysis=isrf_class_analysis
ENDIF
FOR j=0L,Nseds-1 DO BEGIN
;seed=1
inn=abs(randomu(seed))
ii=long(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'
fit_isrf_class_str='_isrfclass'+strtrim(use_isrf_class_analysis,2)
table_name_2=dir+model+'_MuseISRF_JWST_G0_YPAH_YVSG_4Phangs'+fit_isrf_class_str+'.fits'
;print,ii
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,best_sed=best_sed,best_grid_sed=best_grid_sed,all_chi2s=all_chi2s)
;print,params,fact
;IF fact GT 2 THEN BEGIN
; print,ii,fact,chi2
; stop
;ENDIF
facts[j]=fact
chi2s[j]=chi2
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
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
win=0L
!p.multi=[0,1,3]
window,win,xsize=500,ysize=900 & win=win+1
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
window,win,xsize=500,ysize=900 & win=win+1
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]*facts,fitted_params[*,1]*facts,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]*facts,fitted_params[*,2]*facts,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
cleanplot
window,win & win=win+1
diff_pah_perc=abs(true_params[*,1]-fitted_params[*,1])/true_params[*,1]*100.
cgplot,facts,diff_pah_perc,psym=4,xtit='fact',ytit='true-fitted Ypah [%]',/xlog,/ylog
cleanplot
window,win & win=win+1
cgplot,chi2s,diff_pah_perc,psym=4,xtit='chi2',ytit='true-fitted Ypah [%]',/xlog,/ylog
stop
END