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