##  BR: improvements to make on the code:

## take out hard-coding of column numbers, reference by name:
###AJR: Done!
### write to fits instead of .dat.



#Ashley J. Ross v1.2
#Generalized routines to calculate ngal vs. nsys and then weights for stellar density and seeing
#requires healpix.py, basic healpix routines translated to python by Ashley J. Ross
#requires star map, provided is
#allstars17.519.9Healpixall128.dat a map of stars with 17.5 < i < 19.9 at Nside = 128 

from math import *
import numpy
import os
import fitsio
import sys
import re

#angles are assumed to be in degrees, equatorial coordinates assumed throughout

#directory for galaxies and random files
#dir = '/Users/ashleyross/fitsfiles/' 
dir = os.getenv('BOSSLSS_DIR') + "/data/"
#directory for systematic maps
#dirsys = '/Users/ashleyr/weightcode/' 
## this should be passed down from mksamplemaster.py, but for now we'll hard-code it.
dirsys = '../inputfiles/'

#parameters that might need to be changed
ic=3 #we are using the i band
zmin=.43 #minimum redshift
zmax=.7 #maximum redshift
app='-Anderson.dat.fits' #this is how all the mksample files end up on my computer
appr='-Anderson.ran.fits'

def luptm(nmag,bnd):
	b = []
	b.append(1.4*10**-10.)
	b.append(.9*10**-10.)
	b.append(1.2*10**-10.)
	b.append(1.8*10**-10.)
	b.append(7.4*10**-10.)
	return -2.5/log(10.)*(asinh((nmag/10.**9.)/(2.*b[bnd]))+log(b[bnd]))

## BR 5/1/14: change default behavior to use the full randoms file.
#def ranHealp_fits(file,NS,res=128,nrx=1000000):
def ranHealp_fits(file,NS,res=128,nrx=-1):
	from healpix import healpix,radec2thphi,thphi2radec
	pixl = []
	h = healpix()
	np = 12*res**2
	for i in range(0,np):
		pixl.append(0)
	if nrx == -1:
		f = fitsio.read(dir+file+'-'+NS+appr)
	else:
		f = fitsio.read(dir+file+'-'+NS+appr)[0:nrx]
	if len(f) == nrx:
		print 'probably not using all available randoms...'
	nru = len(f)
	
	print 'putting '+str(nru)+' randoms into pixel file'	
	print 'this is the filename: ',dir+file+'-'+NS+appr
	for i in range(0,nru):
		#d = f[1][i]			
		#ra,dec = d['RA'][0],d['DEC'][0]
			#print ra,dec
		ra,dec = f[i]['RA'],f[i]['DEC']
		th,phi = radec2thphi(ra,dec)
		p = int(h.ang2pix_nest(res,th,phi))
		pixl[p] += 1.
		if i/100000 == i/100000.:
			print i
	fo = open('ranHeal_pix'+str(res)+'_'+file+'-'+NS+'.dat','w')
	for i in range(0,len(pixl)):
		if pixl[i] > 0:
			fo.write(str(i)+' '+str(pixl[i])+'\n')
	pixl = []
	fo.close()
	return True



def skyseeair2hpix(res=512):
	from healpix import healpix,radec2thphi,thphi2radec	
	pixl = []
	seepl = []
	airpl = []
	nl = []
	np = 12*res**2.
	h = healpix()
	for i in range(0,np):
		pixl.append(0)
		nl.append(0)
		seepl.append(0)
		airpl.append(0)
	f = open('airseesky19.219.8.csv')
	f.readline()
	for line in f:
		ln = line.split(',')
		sky = float(ln[4])
		see = float(ln[3])
		air = float(ln[2])
		ra = float(ln[0])
		dec = float(ln[1])
		th,phi = radec2thphi(ra,dec)
		p = int(h.ang2pix_nest(res,th,phi))
		pixl[p] += sky
		seepl[p] += see
		airpl[p] += air
		nl[p] += 1.
	f = open('airseesky19.820.4.csv')
	f.readline()
	for line in f:
		ln = line.split(',')
		sky = float(ln[4])
		see = float(ln[3])
		air = float(ln[2])
		ra = float(ln[0])
		dec = float(ln[1])
		th,phi = radec2thphi(ra,dec)
		p = int(h.ang2pix_nest(res,th,phi))
		pixl[p] += sky
		seepl[p] += see
		airpl[p] += air
		nl[p] += 1.
	f = open('airseesky20.420.6.csv')
	f.readline()
	for line in f:
		ln = line.split(',')
		sky = float(ln[4])
		see = float(ln[3])
		air = float(ln[2])
		ra = float(ln[0])
		dec = float(ln[1])
		th,phi = radec2thphi(ra,dec)
		p = int(h.ang2pix_nest(res,th,phi))
		pixl[p] += sky
		seepl[p] += see
		airpl[p] += air
		nl[p] += 1.
	for i in range(0,len(pixl)):
		if nl[i] > 0:
			pixl[i] = pixl[i]/nl[i]
			seepl[i] = seepl[i]/nl[i]
			airpl[i] = airpl[i]/nl[i]


	fw  = 'nm'
	fsky = open('healsky'+fw+str(res)+'.dat','w')
	fsee = open('healsee'+fw+str(res)+'.dat','w')
	fair = open('healair'+fw+str(res)+'.dat','w')
		
	for i in range(0,len(pixl)):
		fsky.write(str(pixl[i])+'\n')
		fsee.write(str(seepl[i])+'\n')
		fair.write(str(airpl[i])+'\n')
				
	fsky.close()
	fsee.close()
	fair.close()
	return True

		
def mkpixl_fits(file,res,min=0,max=1,c=999,wm='cp'):
	from healpix import healpix,radec2thphi
	import fitsio
	np = 12*res*res
	f = fitsio.read(dir+file+app) #written for mksample format, is the double fits something I did?
	if wm == 'stm':
	 	fstw = open(file+'-stweights.dat')
	if wm == 'seem':
		fseew = open(file+'-seeweights.dat') 	
	h = healpix()
	pixl = []
	for i in range(0,np):
		pixl.append(0)
	n = 0
	g = 0
	#assumes file is already masked properly and has equatorial ra,dec coordinates in degrees
	for i in range (0,len(f)):
		k = 0
		if wm == 'stm':
			wstm = float(fstw.readline())
		if wm == 'seem':
			wsee = float(fseew.readline())	
		if c != 999:
			if c == 'ifib':
				t = f[i]['FIBER2FLUX'][ic] 
				fibi = luptm(t,ic)-f[i]['EXTINCTION'][ic] 
				if fibi < min or fibi > max:
					k = 1
		z = f[i]['Z']
		if z < zmin or z > zmax:
			k = 1		
		if k == 0:
			ra,dec = f[i]['RA'],f[i]['DEC']
			w = 1. #if there are weights to be used in the file, put them here
			if wm == 'cp':
				w = f[i]['WEIGHT_NOZ']+f[i]['WEIGHT_CP']-1.
			if wm == 'st':
				w = (f[i]['WEIGHT_NOZ']+f[i]['WEIGHT_CP']-1.)*f[i]['WEIGHT_STAR']
			if wm == 'stm':
				w = (f[i]['WEIGHT_NOZ']+f[i]['WEIGHT_CP']-1.)*wstm
			if wm == 'see':
				w = (f[i]['WEIGHT_NOZ']+f[i]['WEIGHT_CP']-1.)*f[i]['WEIGHT_SEEING']
			if wm == 'seem':
				w = (f[i]['WEIGHT_NOZ']+f[i]['WEIGHT_CP']-1.)*wsee			

			th,phi = radec2thphi(ra,dec)
			p = int(h.ang2pix_nest(res,th,phi))
			pixl[p] += 1.*w
			g += 1.
		else:
			n += 1.
	print g,n			
	return pixl

## BR 5/1/14: change default behavior to use the full randoms file.
def ngvsee_ran(file,maskf='',sysmin=.7,sysmax=2.,wm='cp',bnd=3,bnd2='',imin=0,imax=30.,zmmin=0,zmmax=30.,md='see',nrx=-1,mdf=.183):
	#file is string for mksample fits file without 
	#using seeing values from galaxy/randoms files, no pixelization
	nrx = int(nrx)

	binnbs = []
	binns = []
	nsysbin = 10 #10 bins are being used
	for i in range(0,nsysbin):
		binnbs.append(0)
		binns.append(0)

	nbt = 0
	nt = 0
	bs = 0
	bsr = 0
	n0 = 0
	sysm = float(nsysbin)/(sysmax-sysmin)
	ng0 = 0
	
	#make histogram of galaxy numbers vs. seeing
	f = fitsio.read(dir+file+app) #written for mksample format, is the double fits something I did?
	if wm == 'stm':
	 	fstw = open(file+'-stweights.dat')
	if wm == 'seem':
		fseew = open(file+'-seeweights.dat') 	
	n = 0
	g = 0
	for i in range(0,len(f)):
		k = 0
		if wm == 'stm':
			wstm = float(fstw.readline())
		if wm == 'seem':
			wsee = float(fseew.readline())	
		#print d
		z = f[i]['Z']
		if z < zmin or z > zmax:
			k = 1
			
		if bnd2 != '':
			see = abs(f[i]['PSF_FWHM'][bnd]-f[i]['PSF_FWHM'][bnd2]-mdf)
		else:	
			see = f[i]['PSF_FWHM'][bnd]
		#if see < sysmin or see > sysmax:
		#	k = 1			
		imod = luptm(f[i]['MODELFLUX'][3],3)-f[i]['EXTINCTION'][3]
		if imod < imin or imod > imax:
			k = 1	
		zmod = luptm(f[i]['MODELFLUX'][4],4)-f[i]['EXTINCTION'][4]
		if zmod < zmmin or zmod > zmmax:
			k = 1	
		if k == 0:
			w = 1. #if there are weights to be used in the file, put them here
			if wm == 'cp':
				w = f[i]['WEIGHT_NOZ']+f[i]['WEIGHT_CP']-1.
			if wm == 'st':
				w = (f[i]['WEIGHT_NOZ']+f[i]['WEIGHT_CP']-1.)*f[i]['WEIGHT_STAR']
			if wm == 'stm':
				w = (f[i]['WEIGHT_NOZ']+f[i]['WEIGHT_CP']-1.)*wstm
			if wm == 'see':
				w = (f[i]['WEIGHT_NOZ']+f[i]['WEIGHT_CP']-1.)*f[i]['WEIGHT_SEEING']
			if wm == 'seem':
				w = (f[i]['WEIGHT_NOZ']+f[i]['WEIGHT_CP']-1.)*wsee			
			nbt += w
			bins = int((see-sysmin)*sysm)
			if bins >= 0 and bins < nsysbin:
				binnbs[bins] += w
			else:
				bs += w #count numbers outside of sysmin/sysmax
			
			g += 1.
		else:
			n += 1.
	print g,n			
	#make histogram for randoms vs. seeing
	#assumes seeing is the last column in the random file
	#f.close()
	msee = 0
	if nrx == -1:
		f = fitsio.read(dir+file+appr)
	else:
		f = fitsio.read(dir+file+appr)[0:nrx]


	#f = fitsio.FITS(dir+file+appr)
	#nr = f[1].get_nrows() #note, this is much faster than the alternative...
	#if nr < nrx:
	#	nru = nr
	#else:
	#	nru = nrx
	#	print 'not using all available randoms...'
	nru = len(f)
	if nru == nrx:
		print 'probably not using all available randoms...'
	print 'putting '+str(nru)+' randoms into seeing histogram'	
	#d = f[1]['PSF_FWHM'][0:nru]
	for i in range(0,nru):
		d = f[i]['PSF_FWHM']
		if bnd2 != '':
			#see = abs(d[i][bnd]-d[i][bnd2]-mdf)
			see = abs(d[bnd]-d[bnd2]-mdf)
			#ms = d[i][bnd]-d[i][bnd2]
			ms = d[bnd]-d[bnd2]
		else:	
			#see = d[i][bnd]
			see = d[bnd]
			ms = see		
		if see > sysmin and see < sysmax:
			nt += 1.
			bins = int((see-sysmin)*sysm)
			if bins >= 0 and bins < nsysbin:
				binns[bins] += 1.
			msee += ms	
		else:
			bsr += 1. #count numbers outside of sysmin/sysmax
	msee = msee/nt
	print 'mean (delta) seeing '+str(msee)
	print 'total number of randoms(or pixels)/objects '+str(nt)+'/'+str(nbt)
	print 'number of randoms(or pixels)/objects outside tested range '+str(bsr)+'/'+str(bs)			
	ave = nbt/nt
	print 'average number of objects per random (or pixel) is '+ str(ave)
	#fs = open('n'+file+str(res)+str(wm)+'vnsee.dat','w')
	isw = ''
	if imin != 0:
		isw += 'im'+str(imin)
	if imax != 30.:
		isw += 'ix'+str(imax)
	zsw = ''
	if zmmin != 0:
		zsw += 'zmm'+str(zmmin)
	if imax != 30.:
		isw += 'zmx'+str(zmmax)
	fs = open('n'+file+'ran'+str(wm)+str(bnd)+str(bnd2)+isw+zsw+'v'+md+'.dat','w')
	xl = []
	yl = []
	el = []
	for i in range(0,nsysbin):
		sysv = sysmin + 1./(2.*sysm) + i/sysm
		if binns[i] > 0:
			ns = binnbs[i]/binns[i]/ave
			nse = sqrt(binnbs[i]/(binns[i])**2./(ave)**2.+(binnbs[i]/ave)**2./(binns[i])**3.) #calculate poisson error
		else:
			ns = 1. #write out 1.0 1.0 if no pixels at given value of sys
			nse = 1.		
		fs.write(str(sysv)+' '+str(ns)+' '+str(nse)+'\n')
		xl.append(sysv)
		yl.append(ns)
		el.append(nse)
	fs.close()
	return xl,yl,el


def ngvstar_ifib(file,ifibm,ifibx,maskf='',sysmin=100,sysmax=1300,res=128,wm='cp'):
	#file is string for mksample fits file without 
	#hardcoded columns for relevant quantities, probably a better way to write this...
	ic = 3
	rac = 0
	decc = 1
	ifibc = 16
	extic = 17 #question: should we extinction correct ifib...?
	stl = []
	wstl = []
	errl = []
	fsys = open(dirsys+'allstars17.519.9Healpixall'+str(res)+'.dat').readlines()

	ml = []
	npo = 12*res**2
	for i in range(0,npo):
		ml.append(0)
	mm = 0
	f = open('ranHeal_pix'+str(res)+maskf+'.dat') #created using ranHealp_fits
	
	for line in f:
		ln = line.split()
		pix = int(ln[0])
		ml[pix] = float(ln[1])

	pixl = mkpixl_fits(file,res,min=ifibm,max=ifibx,c='ifib',wm=wm)
	binnbs = []
	binns = []
	nsysbin = 10 #10 bins are being used
	for i in range(0,nsysbin):
		binnbs.append(0)
		binns.append(0)

	nbt = 0
	nt = 0
	bs = 0
	bsr = 0
	n0 = 0
	sysm = float(nsysbin)/(sysmax-sysmin)
	ng0 = 0
	for i in range(0,npo):
		sysv = float(fsys[i])
		if sysv != 0: #the maps are not perfect, entries with 0s shouldn't be used
			nt += ml[i]
			nbt += pixl[i]
			bins = int((sysv-sysmin)*sysm)
			if bins >= 0 and bins < nsysbin:
				binnbs[bins] += pixl[i]
				binns[bins] += ml[i]
			else:
				bs += pixl[i] #count numbers outside of sysmin/sysmax
				bsr +=  ml[i]
		else:
			n0 += ml[i] #count numbers inside bad pixels in sys map
			ng0 += pixl[i]
					
	print 'total number of randoms(or pixels)/objects '+str(nt)+'/'+str(nbt)
	print 'number of randoms(or pixels)/objects where sys = 0 '+str(n0)+'/'+str(ng0)
	print 'number of randoms(or pixels)/objects outside tested range '+str(bsr)+'/'+str(bs)			
	ave = nbt/nt
	print 'average number of objects per random (or pixel) is '+ str(ave)
	fs = open('n'+file+str(res)+'ifib'+str(ifibm)+str(ifibx)+str(wm)+'vnst.dat','w')
	xl = []
	yl = []
	el = []
	for i in range(0,nsysbin):
		sysv = sysmin + 1./(2.*sysm) + i/sysm
		if binns[i] > 0:
			ns = binnbs[i]/binns[i]/ave
			nse = sqrt(binnbs[i]/(binns[i])**2./(ave)**2.+(binnbs[i]/ave)**2./(binns[i])**3.) #calculate poisson error
		else:
			ns = 1. #write out 1.0 1.0 if no pixels at given value of sys
			nse = 1.		
		fs.write(str(sysv)+' '+str(ns)+' '+str(nse)+'\n')
		xl.append(sysv)
		yl.append(ns)
		el.append(nse)

	return xl,yl,el
	

def ifiblinfits(v='cmass-dr11v1'):
	#combines NGC+SGC for five ifib bins used in previous papers
	from optimize import fmin
	rn = ngvstar_ifib(v+'-N',0,20.3,maskf='_'+v+'-N')
	rs = ngvstar_ifib(v+'-S',0,20.3,maskf='_'+v+'-S')
	y0t = []
	e0t = []
	for i in range(0,len(rn[1])):
		yt = (rn[1][i]/rn[2][i]**2.+rs[1][i]/rs[2][i]**2.)/(1./rn[2][i]**2.+1./rs[2][i]**2.)
		et = sqrt(1./(1./rn[2][i]**2.+1./rs[2][i]**2.))
		y0t.append(yt)
		e0t.append(et)
	lf = linfit(rn[0],y0t,e0t)
	inl = numpy.array([1.,0])
	b0,m0 = fmin(lf.chilin,inl)
	print y0t[4],rn[1][4],e0t[4],rn[2][4] #check that answers make sense
	print b0,m0
	#return True
	rn = ngvstar_ifib(v+'-N',20.3,20.6,maskf='_'+v+'-N')
	rs = ngvstar_ifib(v+'-S',20.3,20.6,maskf='_'+v+'-S')
	y1t = []
	e1t = []
	for i in range(0,len(rn[1])):
		yt = (rn[1][i]/rn[2][i]**2.+rs[1][i]/rs[2][i]**2.)/(1./rn[2][i]**2.+1./rs[2][i]**2.)
		et = sqrt(1./(1./rn[2][i]**2.+1./rs[2][i]**2.))
		y1t.append(yt)
		e1t.append(et)
	lf = linfit(rn[0],y1t,e1t)
	inl = numpy.array([1.,0])
	b1,m1 = fmin(lf.chilin,inl)
	print b1,m1
	rn = ngvstar_ifib(v+'-N',20.6,20.9,maskf='_'+v+'-N')
	rs = ngvstar_ifib(v+'-S',20.6,20.9,maskf='_'+v+'-S')
	y2t = []
	e2t = []
	for i in range(0,len(rn[1])):
		yt = (rn[1][i]/rn[2][i]**2.+rs[1][i]/rs[2][i]**2.)/(1./rn[2][i]**2.+1./rs[2][i]**2.)
		et = sqrt(1./(1./rn[2][i]**2.+1./rs[2][i]**2.))
		y2t.append(yt)
		e2t.append(et)
	lf = linfit(rn[0],y2t,e2t)
	inl = numpy.array([1.,0])
	b2,m2 = fmin(lf.chilin,inl)
	print b2,m2
	rn = ngvstar_ifib(v+'-N',20.9,21.2,maskf='_'+v+'-N')
	rs = ngvstar_ifib(v+'-S',20.9,21.2,maskf='_'+v+'-S')
	y3t = []
	e3t = []
	for i in range(0,len(rn[1])):
		yt = (rn[1][i]/rn[2][i]**2.+rs[1][i]/rs[2][i]**2.)/(1./rn[2][i]**2.+1./rs[2][i]**2.)
		et = sqrt(1./(1./rn[2][i]**2.+1./rs[2][i]**2.))
		y3t.append(yt)
		e3t.append(et)
	lf = linfit(rn[0],y3t,e3t)
	inl = numpy.array([1.,0])
	b3,m3 = fmin(lf.chilin,inl)
	print b3,m3
	rn = ngvstar_ifib(v+'-N',21.2,30,maskf='_'+v+'-N')
	rs = ngvstar_ifib(v+'-S',21.2,30,maskf='_'+v+'-S')
	y4t = []
	e4t = []
	for i in range(0,len(rn[1])):
		yt = (rn[1][i]/rn[2][i]**2.+rs[1][i]/rs[2][i]**2.)/(1./rn[2][i]**2.+1./rs[2][i]**2.)
		et = sqrt(1./(1./rn[2][i]**2.+1./rs[2][i]**2.))
		y4t.append(yt)
		e4t.append(et)
	lf = linfit(rn[0],y4t,e4t)
	inl = numpy.array([1.,0])
	b4,m4 = fmin(lf.chilin,inl)
	fo = open('nstlinfits'+v+'128.dat','w')
	fo.write(str(b0)+' '+str(m0)+'\n')
	fo.write(str(b1)+' '+str(m1)+'\n')
	fo.write(str(b2)+' '+str(m2)+'\n')
	fo.write(str(b3)+' '+str(m3)+'\n')
	fo.write(str(b4)+' '+str(m4)+'\n')
	fo.close()
	return True

## old header for this function!  stripped out file names to be more generic.
#def assignstweights(v='cmass-dr11v1',NS='N',res=128,rc=0,dc=1):
def assignstweightsgeneric(fitsfname,stparamfile,outtxtfname=None,res=128,rc=0,dc=1):
	from healpix import healpix,radec2thphi
	import fitsio
	np = 12*res*res
	f = fitsio.read(fitsfname) 
	## make a vector to return
	wst = numpy.zeros(len(f))
	h = healpix()
	lfits = open(stparamfile).readlines()   #'nstlinfits'+v+'128.dat').readlines()
	fsys = open(dirsys+'allstars17.519.9Healpixall'+str(res)+'.dat').readlines()
	if(outtxtfname is not None):
		fo = open(outtxtfname,'w')
	for i in range(0,len(f)):
		t = f[i]['FIBER2FLUX'][ic] 
		fibi = luptm(t,ic)-f[i]['EXTINCTION'][ic] 
		if fibi < 20.45:
			b0 = float(lfits[0].split()[0])
			m0 = float(lfits[0].split()[1])
			b1 = float(lfits[1].split()[0])
			m1 = float(lfits[1].split()[1])
			x0 = 20.15
			x1 = 20.45
		if fibi >= 20.45 and fibi < 20.75:
			b0 = float(lfits[1].split()[0])
			m0 = float(lfits[1].split()[1])
			b1 = float(lfits[2].split()[0])
			m1 = float(lfits[2].split()[1])
			x0 = 20.45
			x1 = 20.75
		if fibi >= 20.75 and fibi < 21.05:
			b0 = float(lfits[2].split()[0])
			m0 = float(lfits[2].split()[1])
			b1 = float(lfits[3].split()[0])
			m1 = float(lfits[3].split()[1])
			x0 = 20.75
			x1 = 21.05
		if fibi >= 21.05:
			b0 = float(lfits[3].split()[0])
			m0 = float(lfits[3].split()[1])
			b1 = float(lfits[4].split()[0])
			m1 = float(lfits[4].split()[1])
			x0 = 21.05
			x1 = 21.35
		mm = (m1-m0)/(x1-x0)
		bm = m0-mm*x0
		mb = (b1-b0)/(x1-x0)
		bb = b0-mb*x0
		mf = mm*fibi+bm
		bf = mb*fibi+bb		
		ra,dec = f[i][rc],f[i][dc]
		th,phi = radec2thphi(ra,dec)
		p = int(h.ang2pix_nest(res,th,phi))
		nst = float(fsys[p])
		wst[i] = 1./(mf*nst+bf)
		if(outtxtfname is not None):
			fo.write(str(wst[i])+'\n')
	if(outtxtfname is not None):
		fo.close()
	## was return True! 
	return wst

def assignstweights(v='cmass-dr11v1',NS='N',res=128,rc=0,dc=1):
  fitsfname = dir+v+'-'+NS+app
  outtxtfname = v+'-'+NS+'-stweights.dat'
  stparamfile = 'nstlinfits'+v+'128.dat'
  assert res == 128
  return assignstweightsgeneric(fitsfname,stparamfile,outtxtfname,res,rc,dc)

def assignseeweights(v='cmass-dr11v1',NS='N'):
	f = fitsio.read(dir+v+'-'+NS+app) 
	fo = open(v+'-'+NS+'-seeweights.dat','w')
	A,B,sig = seew(v,NS)
        ## write out the parameters!
        seeparamfile = 'seefits'+v+'-'+NS+'.dat'
        ofp = open(seeparamfile,'w')
        ofp.write('%.15e %.15e %.15e\n' % (A,B,sig))
        ofp.close()
	for i in range(0,len(f)):
		see = f[i]['PSF_FWHM'][3]
		wsee = 1./(A*(1.-erf((see-B)/sig)))
		fo.write(str(wsee)+'\n')
	fo.close()
	fsys = []
	return True

def assignseeweightsLZ36(v='lowzearly3-6-dr12v5',NS='N'):
	from optimize import fmin
	f = fitsio.read(dir+v+'-'+NS+app) 
	fo = open(v+'-'+NS+'-seeweights.dat','w')
	x,y,e = ngvsee_ran(v+'-N',imin=18.5)
	lf = linfit(x,y,e)
	inl = numpy.array([1.,0])
	b0,m0 = fmin(lf.chilin,inl)
	x,y,e = ngvsee_ran(v+'-N',imin=18.,imax=18.5)
	lf = linfit(x,y,e)
	inl = numpy.array([1.,0])
	b1,m1 = fmin(lf.chilin,inl)
	x,y,e = ngvsee_ran(v+'-N',imin=17.5,imax=18.)
	lf = linfit(x,y,e)
	inl = numpy.array([1.,0])
	b2,m2 = fmin(lf.chilin,inl)
	x,y,e = ngvsee_ran(v+'-N',imax=17.5)
	lf = linfit(x,y,e)
	inl = numpy.array([1.,0])
	b3,m3 = fmin(lf.chilin,inl)
	ml = numpy.array([m3,m2,m1,m0])
	xl = numpy.array([17.25,17.75,18.25,18.75])
	x = (xl-16.)**-.5
	e = numpy.ones((4))
	inl = numpy.array([1.,0])
	lf = linfit(x,ml,e)
	b,m = fmin(lf.chilin,inl)
	
	fp = open('seefits'+v+'-'+NS+'.dat','w')
	fp.write('iband intercept slope\n')
	fp.write(str(17.25)+' '+str(b3)+' '+str(m3)+'\n')
	fp.write(str(17.75)+' '+str(b2)+' '+str(m2)+'\n')
	fp.write(str(18.25)+' '+str(b1)+' '+str(m1)+'\n')
	fp.write(str(18.75)+' '+str(b0)+' '+str(m0)+'\n')
	fp.write('intercept slope fit to slp = (imod-16)^-.5:\n')
	fp.write(str(b)+' '+str(m)+'\n')
	fp.close()
	#b0 = b0/(b0+1.25*m0) #makes the weight 1 at see_i = 1.25
        ## write out the parameters!
	#seeparamfile = 'seefits'+v+'-'+NS+'.dat'
	#ofp = open(seeparamfile,'w')
	#ofp.write('%.15e %.15e \n' % (b0,m0))
	#ofp.close()
	for i in range(0,len(f)):
		see = f[i]['PSF_FWHM'][3]
		imod = luptm(f[i]['MODELFLUX'][3],3)-f[i]['EXTINCTION'][3]
		slp = b + m*(imod-16.)**-.5
		if slp < -2.:
			slp = -2.
		bi = 1.-1.25*slp #makes the weight 1 at see_i = 1.25 
		nsee = bi+slp*see
		if nsee < .2: #don't allow weight to be greater than 5
			nsee = .2
		wsee = 1./nsee
		fo.write(str(wsee)+'\n')
	fo.close()
	fsys = []
	return True

def assignseeweightsfromparamsLZ36(fitsfname,b,m,outtxtfname=None):
	f = fitsio.read(fitsfname) 
	wsee = numpy.zeros(len(f))
        if(outtxtfname is not None):
		fo = open(outtxtfname,'w')

	for i in range(0,len(f)):
		see = f[i]['PSF_FWHM'][3]
		imod = luptm(f[i]['MODELFLUX'][3],3)-f[i]['EXTINCTION'][3]
		slp = b + m*(imod-16.)**-.5
		if slp < -2.:
			slp = -2.
		bi = 1.-1.25*slp #makes the weight 1 at see_i = 1.25 
		nsee = bi+slp*see
		if nsee < .2: #don't allow weight to be greater than 5
			nsee = .2
		wsee[i] = 1./nsee
        	if(outtxtfname is not None):
			fo.write(str(wsee[i])+'\n')
        if(outtxtfname is not None):
		fo.close()
	return wsee


def assignseeweightsfromparams(fitsfname,A,B,sig,outtxtfname=None):
	f = fitsio.read(fitsfname) 
	wsee = numpy.zeros(len(f))
        if(outtxtfname is not None):
		fo = open(outtxtfname,'w')

	for i in range(0,len(f)):
		see = f[i]['PSF_FWHM'][3]
		wsee[i] = 1./(A*(1.-erf((see-B)/sig)))
        	if(outtxtfname is not None):
			fo.write(str(wsee[i])+'\n')
        if(outtxtfname is not None):
		fo.close()
	fsys = []
	return wsee

def assignseeweightsfromparamfile(fitsfname,seeparamfile,sample,outtxtfname=None):
	if sample == 'lowzearly3-6':
		ifp = open(seeparamfile,'r')
		ll = ifp.readlines()
		ifp.close()
		b = float(ll[6].split(' ')[0])
		m = float(ll[6].split(' ')[1])
		print 'seeing weights for LZ36 params',b,m
		return assignseeweightsfromparamsLZ36(fitsfname,b,m,outtxtfname=outtxtfname)
	else:
		assert sample == 'cmass'
		A,B,sig = numpy.loadtxt(seeparamfile,unpack=True)
        	return assignseeweightsfromparams(fitsfname,A,B,sig,outtxtfname=outtxtfname)
	



	
class linfit:
	def __init__(self,xl,yl,el):
		self.xl = xl
		self.yl = yl
		self.el = el
		
	def chilin(self,bml):
		chi = 0
		b = bml[0]
		m = bml[1]
		for i in range(0,len(self.xl)):
			y = b+m*self.xl[i]
			chi += (self.yl[i]-y)**2./self.el[i]**2.
		return chi	

def seew(file,NS):
	from optimize import fmin
	sf = seefit(file,NS)
	inl = numpy.array([1.,2.,.8])
	A,B,sig = fmin(sf.chitest,inl)
	return A,B,sig
	

class seefit:
	def __init__(self,file,NS):
		self.xl = []
		self.yl = []
		self.el = []
		f = open('n'+file+'-'+NS+'rancp3vsee.dat').readlines()
		#fs = open('n'+file+'-S'+'rancp3vsee.dat').readlines()
		for i in range(0,len(f)):
			self.xl.append(float(f[i].split()[0]))
			self.yl.append(float(f[i].split()[1]))
			self.el.append(float(f[i].split()[2]))
		
		
	def chitest(self,bml):
		chi = 0
		A = bml[0]
		B = bml[1]
		sig = bml[2]
		for i in range(0,len(self.xl)):
			y = A*(1.-erf((self.xl[i]-B)/sig))
			chi += (self.yl[i]-y)**2./self.el[i]**2.
		return chi	
			
			
	

def chitestplot_star(file,ifm,ifx,xlab=r'$n_{star}/Nside=128$',wm='stm',sysmin=100,sysmax=1300):
	import matplotlib.pyplot as plt
	from matplotlib import rc
	from matplotlib.backends.backend_pdf import PdfPages
	pp = PdfPages('n'+file+str(ifm)+str(ifx)+'vnstar.pdf')
	plt.clf()
	#rc('text',usetex=True)
	#rc('font',**{'family':'serif','serif':['Computer Modern']})
	#rc('font',**{'family':'serif','serif':['STIX sans-serif']})
	#plt.ylabel(r'$n/\bar n$',size=16)
	title=file+' '+str(ifm)+' < ifib < '+str(ifx)
	plt.minorticks_on()
	plt.ylabel('number density / average number density',size=16)
	plt.xlabel(xlab,size=16)
	xnl,ynl,enl = ngvstar_ifib(file+'-N',ifm,ifx,maskf='_'+file+'-N',sysmin=sysmin,sysmax=sysmax,wm=wm)
	chin = 0
	for i in range(0,len(xnl)):
		chin += (ynl[i]-1.)**2./enl[i]**2.
	print 'chi2 is '+str(chin)+' for 10 dof for NGC for '+str(ifm)+' < ifib < '+str(ifx)
	xsl,ysl,esl = ngvstar_ifib(file+'-S',ifm,ifx,maskf='_'+file+'-S',sysmin=sysmin,sysmax=sysmax,wm=wm)
	chis = 0
	for i in range(0,len(xsl)):
		chis += (ysl[i]-1.)**2./esl[i]**2.
	print 'chi2 is '+str(chis)+' for 10 dof for SGC for '+str(ifm)+' < ifib < '+str(ifx)
	ytl = []
	etl = []
	for i in range(0,len(xnl)):
		yt = (ynl[i]/enl[i]/enl[i]+ysl[i]/esl[i]/esl[i])/(1./enl[i]/enl[i]+1./esl[i]/esl[i])
		ytl.append(yt)
		et = sqrt(1./(1./enl[i]/enl[i]+1./esl[i]/esl[i]))
		etl.append(et)
	chit = 0
	for i in range(0,len(xsl)):
		chit += (ytl[i]-1.)**2./etl[i]**2.
	print 'chi2 is '+str(chit)+' for 10 dof for NGC+SGC for '+str(ifm)+' < ifib < '+str(ifx)	
	ol = []
	for i in range(0,len(xnl)):
		ol.append(1.)
	plt.ylim(.801,1.2)
	plt.xlim(sysmin,sysmax)
	#plt.errorbar(xl, yl,el,fmt=None)
	plt.plot(xnl,ol,'--',color='k')
	## didn't work on riemann: TypeError: There is no line property "capthick"
	#plt.errorbar(xnl,ynl,enl,fmt='bo',markersize=7,elinewidth=2,capthick=2,capsize=4,color='r')
	#plt.errorbar(xnl,ysl,esl,fmt='bo',markersize=7,elinewidth=2,capthick=2,capsize=4,color='b')
	#plt.errorbar(xnl,ytl,etl,fmt='bo',markersize=7,elinewidth=2,capthick=2,capsize=4,color='k')
	plt.errorbar(xnl,ynl,enl,fmt='bo',markersize=7,elinewidth=2,color='r')
	plt.errorbar(xnl,ysl,esl,fmt='bo',markersize=7,elinewidth=2,color='b')
	plt.errorbar(xnl,ytl,etl,fmt='bo',markersize=7,elinewidth=2,color='k')
	#plt.plot(xl, yl,'bo',markersize=5)
	plt.text(250,1.15,'NGC+SGC',fontsize=18)
	plt.text(250,1.12,'NGC',fontsize=18,color='r')
	plt.text(250,1.09,'SGC',fontsize=18,color='b')

	plt.title(title)
	pp.savefig()
	pp.close()
	return True

def chitestplot_see(file,xlab='i-band seeing',wm='seem',sysmin=.7,sysmax=2.):
	import matplotlib.pyplot as plt
	from matplotlib import rc
	from matplotlib.backends.backend_pdf import PdfPages
	pp = PdfPages('n'+file+'vsee.pdf')
	plt.clf()
	#rc('text',usetex=True)
	#rc('font',**{'family':'serif','serif':['Computer Modern']})
	#rc('font',**{'family':'serif','serif':['STIX sans-serif']})
	#plt.ylabel(r'$n/\bar n$',size=16)
	title=file
	plt.minorticks_on()
	plt.ylabel('number density / average number density',size=16)
	plt.xlabel(xlab,size=16)
	xnl,ynl,enl = ngvsee(file+'-N',maskf='_'+file+'-N',sysmin=sysmin,sysmax=sysmax,wm=wm)
	chin = 0
	for i in range(0,len(xnl)):
		chin += (ynl[i]-1.)**2./enl[i]**2.
	print 'chi2 is '+str(chin)+' for 10 dof for NGC'
	xsl,ysl,esl = ngvsee(file+'-S',maskf='_'+file+'-S',sysmin=sysmin,sysmax=sysmax,wm=wm)
	chis = 0
	for i in range(0,len(xsl)):
		chis += (ysl[i]-1.)**2./esl[i]**2.
	print 'chi2 is '+str(chis)+' for 10 dof for SGC'
	ytl = []
	etl = []
	for i in range(0,len(xnl)):
		yt = (ynl[i]/enl[i]/enl[i]+ysl[i]/esl[i]/esl[i])/(1./enl[i]/enl[i]+1./esl[i]/esl[i])
		ytl.append(yt)
		et = sqrt(1./(1./enl[i]/enl[i]+1./esl[i]/esl[i]))
		etl.append(et)
	chit = 0
	for i in range(0,len(xsl)):
		chit += (ytl[i]-1.)**2./etl[i]**2.
	print 'chi2 is '+str(chit)+' for 10 dof for NGC+SGC'	
	ol = []
	for i in range(0,len(xnl)):
		ol.append(1.)
	plt.ylim(.801,1.2)
	plt.xlim(sysmin,sysmax)
	#plt.errorbar(xl, yl,el,fmt=None)
	plt.plot(xnl,ol,'--',color='k')
	#plt.errorbar(xnl,ynl,enl,fmt='bo',markersize=7,elinewidth=2,capthick=2,capsize=4,color='r')
	#plt.errorbar(xnl,ysl,esl,fmt='bo',markersize=7,elinewidth=2,capthick=2,capsize=4,color='b')
	#plt.errorbar(xnl,ytl,etl,fmt='bo',markersize=7,elinewidth=2,capthick=2,capsize=4,color='k')
	plt.errorbar(xnl,ynl,enl,fmt='bo',markersize=7,elinewidth=2,color='r')
	plt.errorbar(xnl,ysl,esl,fmt='bo',markersize=7,elinewidth=2,color='b')
	plt.errorbar(xnl,ytl,etl,fmt='bo',markersize=7,elinewidth=2,color='k')
	#plt.plot(xl, yl,'bo',markersize=5)
	plt.text(.9,1.15,'NGC+SGC',fontsize=18)
	plt.text(.9,1.12,'NGC',fontsize=18,color='r')
	plt.text(.9,1.09,'SGC',fontsize=18,color='b')

	plt.title(title)
	pp.savefig()
	pp.close()
	return True


def fillfits(v='cmass-dr11v1',NS='N',input_app='-Anderson.dat.fits',sample='cmass'):
	import fitsio

	global app, appr #zmin,zmax
	## reassign global variables.
        app = input_app

	if sample == 'lowzearly3-6':
        	seewgt = numpy.loadtxt(v+'-'+NS+'-seeweights.dat',unpack=True)
		stwgt = seewgt.copy()
		stwgt[:] = 1.0
	else:
	        stwgt = numpy.loadtxt(v+'-'+NS+'-stweights.dat',unpack=True)
        	seewgt = numpy.loadtxt(v+'-'+NS+'-seeweights.dat',unpack=True)
	print 'opening this file to write systematic weights',dir+v+'-'+NS+app
        dd = fitsio.read(dir+v+'-'+NS+app)
        assert len(dd) == len(stwgt)
        assert len(dd) == len(seewgt)

        dd['WEIGHT_STAR'][:] = stwgt
        dd['WEIGHT_SEEING'][:] = seewgt
        dd['WEIGHT_SYSTOT'][:] = dd['WEIGHT_STAR'][:] * dd['WEIGHT_SEEING'][:]

        print 'about to write to',dir+v+'-'+NS+app
  
        fitsio.write(dir+v+'-'+NS+app,dd,clobber=True)  

def fillfitsfromparams(fitsfname,stparamfile,seeparamfile,sample,writeopt=1):

  seewgt = assignseeweightsfromparamfile(fitsfname,seeparamfile,sample)
  dd = fitsio.read(fitsfname)
  assert len(dd) == len(seewgt)

  if sample == 'lowzearly3-6':
    dd['WEIGHT_STAR'][:] = 1.0

  else:  # standard cmass 
    stwgt = assignstweightsgeneric(fitsfname,stparamfile)
    assert len(dd) == len(stwgt)
    dd['WEIGHT_STAR'][:] = stwgt

  dd['WEIGHT_SEEING'][:] = seewgt
  dd['WEIGHT_SYSTOT'][:] = dd['WEIGHT_STAR'][:] * dd['WEIGHT_SEEING'][:]

  if writeopt == 1:
    fitsio.write(fitsfname,dd,clobber=True)
  return dd


def photokey(darr):
  """
  SDSS magic photoid is some combo of these things, so they can be used as a unique key for each object.
  """
  return (darr['RUN'], darr['RERUN'], darr['CAMCOL'], darr['FIELD'], darr['ID'])


def fillbyObjID(fitsfiletofill,fitsfilefillfrom):
  """
  Nevermind, I don't think this will work in all cases because of the random subsampling of KNOWN objects.
  """
  dtmp = fitsio.read(fitsfilefillfrom)
  targ_objid = {}
  for i in range(len(dtmp)):
    ## see type trickiness in RERUN in match_sample.py
    targ_objid[(dtmp['RUN'][i], dtmp['RERUN'][i], dtmp['CAMCOL'][i], dtmp['FIELD'][i], dtmp['ID'][i])] = i
  dd = fitsio.read(fitsfiletofill)
  for i in range(len(dd)):
    photoid = photokey(dd)
    if photoid in targ_objid:
      dtmpi = targ_objid[photoid]
      dd['WEIGHT_STAR'][i] = dtmp['WEIGHT_STAR'][dtmpi]
      dd['WEIGHT_SEEING'][i] = dtmp['WEIGHT_SEEING'][dtmpi]
      dd['WEIGHT_SYSTOT'][i] = dtmp['WEIGHT_SYSTOT'][dtmpi]
    else:
      print 'bad run/rerun/camcol/field/id match!  Aborting.'
      print dd['RUN'][i], dd['RERUN'][i], dd['CAMCOL'][i], dd['FIELD'][i], dd['ID'][i]
      sys.exit(1)
  del dtmp
  fitsio.write(fitsfiletofill,dd,clobber=True) 

def mkweightfilesLZ36(file,ranopt=1,plotopt=0,input_app='-Anderson.dat.fits',input_appr='-Anderson.ran.fits',input_zmin=0.2,input_zmax=0.5):
	"""
	ranopt = 1 if you have not made healpix random files yet, 0 otherwise
	plotopt = 1 if you want to make plots only, probably need to run this in non-batch mode.
	plotopt = 1 assumes youve already generated the weights files.
	"""
	import sys
	from healpix import ranpixl_up

	global app,appr,zmin,zmax

	## reassign global variables.
	app = input_app
	appr = input_appr
	zmin = input_zmin
	zmax = input_zmax

	
	#find galaxy/seeing relationship for NGC (chunks 3-6 are only in NGC)
	dseefit = 0 #make non-zero if you have done this already
	if dseefit == 0 and plotopt == 0:
		#x,y,e = ngvsee_ran(file+'-N')
		assignseeweightsLZ36(file)
	#write out seeing weights for NGC
	fillfits(file,NS='N',input_app=input_app,sample='lowzearly3-6')

def see36test(file,input_app='-Anderson.dat.fits',input_appr='-Anderson.ran.fits',input_zmin=0.2,input_zmax=0.5):
	import sys
	from healpix import ranpixl_up
	from optimize import fmin
	global app,appr,zmin,zmax

	## reassign global variables.
	app = input_app
	appr = input_appr
	zmin = input_zmin
	zmax = input_zmax
	x,y,e = ngvsee_ran(file+'-N',imin=18.5)
	lf = linfit(x,y,e)
	inl = numpy.array([1.,0])
	b0,m0 = fmin(lf.chilin,inl)
	x,y,e = ngvsee_ran(file+'-N',imin=18.,imax=18.5)
	lf = linfit(x,y,e)
	inl = numpy.array([1.,0])
	b1,m1 = fmin(lf.chilin,inl)
	x,y,e = ngvsee_ran(file+'-N',imin=17.5,imax=18.)
	lf = linfit(x,y,e)
	inl = numpy.array([1.,0])
	b2,m2 = fmin(lf.chilin,inl)
	x,y,e = ngvsee_ran(file+'-N',imax=17.5)
	lf = linfit(x,y,e)
	inl = numpy.array([1.,0])
	b3,m3 = fmin(lf.chilin,inl)
	fo = open('see36mb.dat','w')
	fo.write(str(17.25)+' '+str(b3)+' '+str(m3)+'\n')
	fo.write(str(17.75)+' '+str(b2)+' '+str(m2)+'\n')
	fo.write(str(18.25)+' '+str(b1)+' '+str(m1)+'\n')
	fo.write(str(18.75)+' '+str(b0)+' '+str(m0)+'\n')
	fo.close()
	return True

def test36weight(ep=-2):
	from optimize import fmin
	import numpy
	f = numpy.loadtxt('see36mb.dat').transpose()
	inl = numpy.array([1.,0])
	e = numpy.ones((4))
	x = (f[0]-16.)**ep
	y = f[2]
	inl = numpy.array([1.,0])
	lf = linfit(x,y,e)
	m,b = fmin(lf.chilin,inl)
	return m,b


def mkweightfiles(file,ranopt=1,plotopt=0,input_app='-Anderson.dat.fits',input_appr='-Anderson.ran.fits',input_zmin=0.43,input_zmax=0.7):
	"""
	ranopt = 1 if you have not made healpix random files yet, 0 otherwise
	plotopt = 1 if you want to make plots only, probably need to run this in non-batch mode.
	plotopt = 1 assumes youve already generated the weights files.
	"""
	import sys
	from healpix import ranpixl_up

	global app,appr,zmin,zmax

	## reassign global variables.
	app = input_app
	appr = input_appr
	zmin = input_zmin
	zmax = input_zmax

	#f = fitsio.FITS(dir+file+'-'+NS+appr)
 	#fnameN = dir + filebase + '-N' + app
 	#fnameS = dir + filebase + '-S' + app
	## want to change the code to pass the data vector along, but that is more long term.  first get it working.
	#ddN = fitsio.read(fnameN)
	#ddS = fitsio.read(fnameS)

	#make random pixel files
	#maximum number of randoms set to one million, no reason not use all except for time 
#	ran = 0 #make non-zero if you have already made the random files
## BR: this section creates ranHeal_pix1024* and ranHeal_pix128* files.
## Leave this unmodified.
	if ranopt == 1:
		### BR 5/1/14: changed default behavior to use all the randoms in the input catalog.
		#ranHealp_fits(file,'N',nrx=int(1e6)) 
		#ranHealp_fits(file,'S',nrx=int(1e6)) 
		ranHealp_fits(file,'N') 
		ranHealp_fits(file,'S') 
	
	#find linear fits to ifib binned galaxy/star relationships
	stfits = 0 #make non-zero if you have done these already
## BR: writes out file like nstlinfitscmass-dr11v1remake128.dat. 
	if stfits == 0 and plotopt == 0:
		ifiblinfits(file)
	
	#write out stellar weights for NGC and SGC
	stweights = 0 #make non-zero if you have made these already
## BR:  writes stweights to a dat file instead of overwriting the fits file. fo = open(v+'-'+NS+'-stweights.dat','w')
	if stweights == 0 and plotopt == 0:
		assignstweights(file,'N')
		assignstweights(file,'S')
	
	#find galaxy/seeing relationship for NGC and then SGC
	dseefit = 0 #make non-zero if you have done this already
	if dseefit == 0 and plotopt == 0:
		ngvsee_ran(file+'-N')
		ngvsee_ran(file+'-S')
	
	#write out seeing weights for NGC and then SGC
	seeweights = 0 #make non-zero if you have made these already
## BR: writes to fo = open(v+'-'+NS+'-seeweights.dat','w')
	if seeweights == 0 and plotopt == 0:
		assignseeweights(file)
		assignseeweights(file,'S')

        if plotopt == 0:
		fillfits(file,NS='N',input_app=input_app)
		fillfits(file,NS='S',input_app=input_app)
		
	if plotopt == 1:
		testst = 0 #make non-zero if you do not want to test stellar density results
		if testst == 0:
			#set values to split on ifib for testing
			ifb1 = 20.5
			ifb2 = 21.
			#make plots, print chi2; no cosmic variance so large samples will have large chi2
			chitestplot_star(file,0,ifb1)
			chitestplot_star(file,ifb1,ifb2)
			chitestplot_star(file,ifb2,30)	
		testsee = 0 #make non-zero if you do not want to test seeing weights
		if testsee == 0:	
			chitestplot_see(file)
		
	
	
	



if __name__ == '__main__':
	import sys
	from healpix import ranpixl_up
	if len(sys.argv) != 4:
		print 'Usage: python mkweightfiles.py filebase ranopt plotopt'
		print 'ranopt = 1 if you have not made healpix random files yet, 0 otherwise'
		print 'plotopt = 1 if you want to make plots only, probably need to run this in non-batch mode.'
		print 'plotopt = 1 assumes youve already generated the weights files.'
		sys.exit(1)

	file = str(sys.argv[1]) #for instance, 'cmass-dr11v1'
        ranopt = int(sys.argv[2])
        plotopt = int(sys.argv[3])

	## moved everything that was in main to a function that I can call from other python code.
        mkweightfiles(file=file,ranopt=ranopt,plotopt=plotopt)
	## yes, input variables get passed to mkweightfiles and globals seem to be set, should work!
        #mkweightfiles(file=file,ranopt=ranopt,plotopt=plotopt,input_app='BAR.dat.fits',input_appr='BAG.ran.fits',input_zmin=-39,input_zmax=892.)