PRO  systematics_weight,match_file,output_file,params_star_file,params_seeing_file

;------------------------------------------------------------------------------
;+
; NAME:
;  systematics_weight
;
; PURPOSE:
;   Calculate the w_sys (systematic weights) which is basically a stellar weight 
;   according to Anderson et al. and Ross et al. *(but mostly Ariel's private communications) 
;   for large-scale structure analysis (DR9).
;   DR10 update includes seeing and stellar density weights.
;
; CALLING SEQUENCE:
;   stellar_weight,input-file (for example collided-cmass.fits) ,output_file 
;
; OPTIONAL INPUTS:
;
; OUTPUTS:
;   collided-cmass-star_weight.fits
;   -- in the ../data directory.
;                                
; COMMENTS:
;
; EXAMPLES:
;
; REQUIREMENTS:
; 
; REVISION HISTORY:
;   20-Apr-2012 Written by Shirley Ho (LBL/CMU)
;   26-sep-2012 Edited by Shirley Ho (CMU)
;   29-Apr-2013 Edited by Beth Reid (LBL)
;-

print,' '
print,'Starting weighting with both stars and seeing.'
print,'Match  file is ',match_file
print,'Output file is ',output_file
print,'params_star file is ',params_star_file
print,'params_seeing file is ',params_seeing_file
print,' '


; read in the constants from the file provided by Ariel: 
;;; read in the file from the command line so we dont mix up dr10/dr11vX versions of the paramaters.
;star_weight_file = 'inputfiles/params_wstar_F_128_dr11v2.dat'
;seeing_weight_file = 'inputfiles/params_wseeing_F_1024_dr11v2.dat'

star_weight_file = params_star_file
seeing_weight_file = params_seeing_file

;; these are more stable, leave them hard-coded for now.
starsmapf = 'inputfiles/nstar_healpix_nside128.dat'
seeingmapf = 'inputfiles/healseenm1024.dat'
load_star_coef, star_weight_file, weightstar,n1

ifib2_grid = weightstar[*,2] 
Astar_grid = weightstar[*,3] 
Bstar_grid = weightstar[*,5]

load_seeing_coef, seeing_weight_file, weightseeing,n2
if(n2 NE 1) then begin
  print,'n2 should be 1!',n2
  stop
endif
Aseeing = weightseeing[0,0]
Bseeing = weightseeing[0,2]
Cseeing = weightseeing[0,4]
print,'seing params:',Aseeing,Bseeing,Cseeing

;input
data_dir  = filepath(root=getenv("BOSSLSS_DIR"),"/data/")
; We can calculate this at any time during the process, since the input i_fib and the stellar densities do not change according to cuts in the data. I just choose to do this after matching the galaxies to the redshfits.
galaxyf = data_dir+match_file

resstar =7L ; nside=128: that's what Ariel used to get the coefficients
resseeing =10L ; nside=1024: that's what Ariel used to get the coefficients
;input systematic map
;starsf = 'nstar_healpix_m19.9_nside128_ascii.dat'  ; number of stars per square degree
;seeingf = 'seeing_healpix_m19.9_nside128_ascii.dat'  ; seeing in arcsec
; I switched to Ariel's map, since this is the map that was used to calibrate the stellar weights


;output
galfits = data_dir+output_file

;reading in the galaxy and systematic maps
gal = mrdfits(galaxyf,1) 
read_systematicsmap,starsmapf,stars,nstars
print,'There are ',nstars,' stars pixels'
read_systematicsmap,seeingmapf,seeing,nseeing
print,'There are ',nseeing,' seeing pixels'

;npix =healnpix(res)
;nside = healnpix(res,/nside)
nsidestar = 2L^resstar
npixstar = 12L*nsidestar*nsidestar
area_pixstar = 4*!PI*(180/!pi)*(180/!PI)/npixstar
print,'area of each pix', area_pixstar, 'square degree'
area_pixstar = area_pixstar*60*60
print,'area of each pix is', area_pixstar, 'square arcminute'

;; copy for seeing
nsideseeing = 2L^resseeing
npixseeing = 12L*nsideseeing*nsideseeing
area_pixseeing = 4*!PI*(180/!pi)*(180/!PI)/npixseeing
print,'area of each pix', area_pixseeing, 'square degree'
area_pixseeing = area_pixseeing*60*60
print,'area of each pix is', area_pixseeing, 'square arcminute'

print,'Turning RA, DEC of each galaxies into ipix -> pixel number'

RA = gal.ra  ; in degrees
DEC = gal.dec 
theta = 90L-DEC
phi = RA
theta =theta/180.*!Pi
phi = phi/180.*!Pi
ang2pix_nest,nsidestar,theta,phi,indexstar ; now that index contains the healpix pixel number (in nested ordering) of each galaxy

;; do same thing but for seeing map.
;; new territory for beth -- if i get the wrong answer, ask Shirley/Martin for help here.
nside = 2L^resseeing
ang2pix_nest,nsideseeing,theta,phi,indexseeing ; now that index contains the healpix pixel number (in nested ordering) of each galaxy

ngal=n_elements(gal.ra)
fiber2mag=fltarr(ngal,5)
for i=0,4 do begin
 fiber2mag[*,i]= 22.5-2.5*alog10(gal.fiber2flux[i]>0.001)-gal.extinction[i]
endfor

Astar_interpol = interpol(Astar_grid,ifib2_grid,fiber2mag[*,3])
Bstar_interpol = interpol(Bstar_grid,ifib2_grid,fiber2mag[*,3])

;; temporarily output the weights I compute to compare against ariel.
;; they seem to match!
;close,11
;tmpfname = 'chkweights.dat'
;openw,11,tmpfname

for i=0L,ngal-1L do begin 
   indstar = indexstar[i]
   indseeing = indexseeing[i]
;we want the stellar weight of each galaxy calculated
   starden = stars[indstar]  ; the star density
   seeingval = seeing[indseeing]  ; the seeing value
   ;; From Ariel's email::
   ;; That is, if a given pixel has seeing > 2.5, the corresponding weight
   ;;w_seeing is computed using S=2.5.
   ;;
   if(seeingval GT 2.5) then begin
     seeingval = 2.5d0
   endif
   gal[i].weight_star = Astar_interpol[i] + (Bstar_interpol[i]*starden)
   gal[i].weight_seeing = 2.0d0/Aseeing/(1.-erf((seeingval-Bseeing)/Cseeing))
   gal[i].weight_systot = gal[i].weight_star*gal[i].weight_seeing
;   printf,11,format='(F,F,F,E,E,E)',gal[i].ra, gal[i].dec, seeingval, gal[i].weight_star, gal[i].weight_seeing, gal[i].weight_systot
endfor   
;close,11
mwrfits,gal,galfits,/create


END