Some IDL code for Meszaros et al. (2015)

README.rst

This gist contains some code used in Meszaros et al. (2015) pertaining to the analysis of multiple populations in globular clusters using the extreme deconvolution algorithm. The necessary data files are available upon request.

The basic use of the code is to run:

IDL> analyze_gc, gc=2, ngauss=2, outfilename='m2pho_2comp.fits'

This runs the analysis on a given globular cluster (see read_gc.pro for the mapping between gc and globular cluster). The produced output file is a fits file that has the ID, the component that the XD code assigns (xdpop), the probability of that assignment (logp, actually the natural log of the probability), and the population that was in the original data file from a preliminary analysis of membership probability. Information on the best-fit amplitudes, means, and covariances of the populations is printed directly.

analyze_gc.pro

PRO ANALYZE_GC, gc=gc, almg=almg, ngauss=ngauss, outfilename=outfilename
;; gc= 2, 3, 5, 13, 15, 92
;; almg= if set, only fit Al Mg
;; ngauss= # of Gaussians
IF ~keyword_set(ngauss) THEN ngauss= 3
;;Read
read_gc, id, abu, pop , gc=gc;;abu = [Mg, Al, Si, Ca, Ti]
IF keyword_set(almg) THEN abu= abu[0:1,*]

;;Prepare for XD
ydata= abu
ycovar= abu*0. ;;errors are uncorrelated, so only need to specify diagonal
IF gc EQ 2 THEN BEGIN
    ycovar[0,*]= 0.08^2.
    ycovar[1,*]= 0.11^2.
    IF ~keyword_set(almg) THEN BEGIN
        ycovar[2,*]= 0.09^2.
        ycovar[3,*]= 0.06^2.
        ycovar[4,*]= 0.13^2.
    ENDIF
ENDIF ELSE IF gc EQ 3 THEN BEGIN
    ycovar[0,*]= 0.04^2.
    ycovar[1,*]= 0.11^2.
    IF ~keyword_set(almg) THEN BEGIN
        ycovar[2,*]= 0.09^2.
        ycovar[3,*]= 0.06^2.
        ycovar[4,*]= 0.16^2.
    ENDIF
ENDIF ELSE IF gc EQ 5 THEN BEGIN
    ycovar[0,*]= 0.05^2.
    ycovar[1,*]= 0.08^2.
    IF ~keyword_set(almg) THEN BEGIN
        ycovar[2,*]= 0.08^2.
        ycovar[3,*]= 0.06^2.
        ycovar[4,*]= 0.13^2.
    ENDIF
ENDIF ELSE IF gc EQ 13 THEN BEGIN
    ycovar[0,*]= 0.06^2.
    ycovar[1,*]= 0.13^2.
    IF ~keyword_set(almg) THEN BEGIN
        ycovar[2,*]= 0.08^2.
        ycovar[3,*]= 0.07^2.
        ycovar[4,*]= 0.16^2.
    ENDIF
ENDIF ELSE IF gc EQ 15 THEN BEGIN
    ycovar[0,*]= 0.09^2.
    ycovar[1,*]= 0.17^2.
    IF ~keyword_set(almg) THEN BEGIN
        ycovar[2,*]= 0.13^2.
        ycovar[3,*]= 0.19^2.
        ycovar[4,*]= 0.28^2.
    ENDIF
ENDIF ELSE IF gc EQ 92 THEN BEGIN
    ycovar[0,*]= 0.09^2.
    ycovar[1,*]= 0.13^2.
    IF ~keyword_set(almg) THEN BEGIN
        ycovar[2,*]= 0.12^2.
        ycovar[3,*]= 0.22^2.
        ycovar[4,*]= 0.20^2.
    ENDIF
ENDIF
;;find bad values and give them large errors to make them essentially
;;unobserved
indx= where(ydata gt 9998.0,cnt)
if cnt gt 0 then begin
    ydata[indx]= 0.
    ycovar[indx]= 100.^2.
endif
IF keyword_set(almg) THEN dx= 2 ELSE dx= 5;;number of abundances

;;Initialize XD
fit_amp= dblarr(ngauss)
fit_mean= dblarr(dx,ngauss)
fit_covar= dblarr(dx,dx,ngauss)
datamean= MEAN(ydata,dimension=2)
datastd= STDDEV(ydata,dimension=2)
FOR ii=0L, ngauss-1 DO BEGIN
    fit_amp[ii]= 1./double(ngauss)
    fit_mean[*,ii]= randomn(seed,dx)*datastd+datamean
ENDFOR
FOR ii= 0L, dx-1 DO fit_covar[ii,ii,*]= 4.*datastd[ii]^2.

;;Run XD
projected_gauss_mixtures_c, ngauss, ydata, ycovar, $
      fit_amp, fit_mean, fit_covar, /quiet
;;print results
print, FORMAT='("Amplitudes of the ",I0," components fit with XD")', ngauss
print, fit_amp
print, FORMAT='("Means of the ",I0," components: Mg, Al, Si, Ca, Ti")', ngauss
print, fit_mean

;;Calculate membership probabilities
calc_membership_prob, logpost, ydata, ycovar, fit_mean, fit_covar, fit_amp
;;logpost is a [ngauss,ndata] array of ln of relative probabilities
maxp= max(logpost,xdpop,dimension=1)
xdpop= xdpop MOD ngauss ;;xdpop has index of highest probability cluster for each star

;;check that we get similar cluster assignments as Szabolcs
FOR ii=0L, ngauss-1 DO BEGIN
    print, FORMAT= '("Cluster ",I0,": old assignments and probability of assignment:")', ii+1
    indx= where(xdpop eq ii,cnt)
    if cnt gt 0 then begin
        print, pop[indx]
        print, exp(maxp[indx])
    endif
ENDFOR

out= {id:0L,oldpop:0L,xdpop:0L,logp:0D}
out= REPLICATE(out,n_elements(id))
out.id= id
out.oldpop= pop
out.xdpop= xdpop
out.logp= maxp
mwrfits, out, outfilename, /create
END

read_gc.pro

PRO READ_GC, id, abu, pop, gc=gc
IF gc EQ 2 THEN BEGIN
    filename= '../data/m2pho.dat'
ENDIF ELSE IF gc EQ 3 THEN BEGIN
    filename= '../data/m3pho.dat'
ENDIF ELSE IF gc EQ 5 THEN BEGIN
    filename= '../data/m5pho.dat'
ENDIF ELSE IF gc EQ 13 THEN BEGIN
    filename= '../data/m13pho.dat'
ENDIF ELSE IF gc EQ 15 THEN BEGIN
    filename= '../data/m15pho.dat'
ENDIF ELSE IF gc EQ 92 THEN BEGIN
    filename= '../data/m92pho.dat'
ENDIF
readcol,filename, format='A,X,D,D,D,D,D,D,D,D,D,D,D,D,L', $
  id, vrad, teff, logg, feh, cfe, nfe, ofe, $
  mgfe, alfe, sife, cafe, tife, pop, skipline=1, delimiter=' '
abu= dblarr(5,n_elements(vrad))
abu[0,*]= mgfe
abu[1,*]= alfe
abu[2,*]= sife
abu[3,*]= cafe
abu[4,*]= tife
END