#!/usr/bin/env python
#Given N-maps calculates cross-reconstruction between them

from flipper import *
from lensTcorrCombIR import *
from numpy.fft import fftshift,fftfreq,fft2,ifft2
from scipy import interpolate
from scipy import *
import os
import random
import sys
import pickle

TCMB = 2.726e6

def getEmptyMapWithDifferentDims(m,Ny,Nx):
   """
   Creates an empty map on the same patch of the sky as m
   but with different dimensions Ny ,Nx

   """
   data = numpy.zeros([Ny,Nx])
   m = m.copy()
   m.wcs.header.update('NAXIS1',Nx)
   m.wcs.header.update('NAXIS2',Ny)
   m.wcs.header.update('CDELT1',  m.wcs.header['CDELT1']*(m.Nx/(Nx*1.0)))
   m.wcs.header.update('CDELT2',  m.wcs.header['CDELT2']*(m.Ny/(Ny*1.0)))
   m.wcs.updateFromHeader()
   p_x, p_y = m.skyToPix(m.x0,m.y0)
   m.wcs.header.update('CRPIX1', m.wcs.header['CRPIX1'] - p_x)
   m.wcs.header.update('CRPIX2', m.wcs.header['CRPIX2'] - p_y)
   m.wcs.updateFromHeader()
   mNew = liteMap.liteMapFromDataAndWCS(data, m.wcs)
   return mNew

def powerFromKMap(k1,k2,ft,ftret):
    """
    @brief Creates a power2D object from one power 2D object and two kmaps
    """
    p2d = ftret.copy()
    mapFFT = k1.copy()
    mapFFT2 = k2.copy()

    area =ft.Nx*ft.Ny*ft.pixScaleX*ft.pixScaleY
    p2d.powerMap = numpy.real(numpy.conjugate(mapFFT)*mapFFT2)*area/(ft.Nx*ft.Ny*1.0)**2
    
    return p2d

def meanCrossSpec(mapList,applySlepianTaper=True,nresForSlepian=3.0):
    count = 0 
    
    for i in xrange(len(mapList)):
        for j in xrange(i):
            
            p2d = fftTools.powerFromLiteMap(mapList[i],mapList[j],\
                                            applySlepianTaper=applySlepianTaper,\
                                            nresForSlepian=nresForSlepian)
            if count == 0:
                p2d0 = p2d.copy()
            else:
                p2d0.powerMap[:] += p2d.powerMap[:]
            count += 1
            
    p2d0.powerMap[:] /= count
    powerM = p2d0.powerMap.copy()
    print 'count=', count
        
    lL,lU,lBin,clBinCrossMean,err,w = p2d0.binInAnnuli(os.environ['LENSRECONS_DIR']+os.path.sep+'params/BIN_100_LOG')
            
    return lBin,clBinCrossMean,powerM

def meanAutoSpec(mapList,applySlepianTaper=True,nresForSlepian=3.0):
    count = 0 
    
    for i in xrange(len(mapList)):
            
        p2d = fftTools.powerFromLiteMap(mapList[i],\
                                            applySlepianTaper=applySlepianTaper,\
                                            nresForSlepian=nresForSlepian)
        if count == 0:
            p2d0 = p2d.copy()
        else:
            p2d0.powerMap[:] += p2d.powerMap[:]
        count += 1
            
    p2d0.powerMap[:] /= count
    powerM2 = p2d0.powerMap.copy()
        
        
    lL,lU,lBin,clBinAutoMean,err,w = p2d0.binInAnnuli(os.environ['LENSRECONS_DIR']+os.path.sep+'params/BIN_100_LOG')
            
    return lBin,clBinAutoMean,powerM2

def meanCrossPowerMap(mapList,applySlepianTaper=True,nresForSlepian=3.0):
    count = 0 
    
    for i in xrange(len(mapList)):
        for j in xrange(i):
            
            p2d = fftTools.powerFromLiteMap(mapList[i],mapList[j],\
                                            applySlepianTaper=applySlepianTaper,\
                                            nresForSlepian=nresForSlepian)
            if count == 0:
                p2d0 = p2d.copy()
            else:
                p2d0.powerMap[:] += p2d.powerMap[:]
            count += 1
            
    p2d0.powerMap[:] /= count

        
    return p2d0


def deconvolveMap(conmap,fwhmRad):

    #deconvolving####################we are only using ls up to 20000 to prevent leakage
    ell = numpy.arange(20000)
    Filter = numpy.exp(ell**2.0*fwhmRad**2.0/(16.*numpy.log(2.)))
    print "deconvolving.."
    ft2 = fftTools.fftFromLiteMap(conmap)
    kMap = ft2.kMap.copy()
    kFilter0 = numpy.real(kMap.copy())*0.+ 1.
    kFilter = kMap.copy()*0.
    FlSpline = splrep(ell,Filter,k=3)
    ll = numpy.ravel(ft2.modLMap)
    kk = (splev(ll,FlSpline))
    kFilter = numpy.reshape(kk,[conmap.Ny,conmap.Nx])
    kFilter0 *= kFilter
    kMap[:,:] *= kFilter0[:,:]
    ft2.kMap = kMap
    deconlensedMap2data = numpy.real(ifft2(kMap))
    conmap.data[:] = deconlensedMap2data[:]
    return conmap

def deconvolveMapRealBeam(freq,ftthing):
    
    if freq == 220:
        ellist = []
        filtlist = []
        files = open('/u/sudeep/ACT/mapAnalysis/compactSources/data/B_l_AR1.dat','r')
        for line in files:
            aa = line.split(' ')
            bb = float(aa[0])
            cc = 1./float(aa[1])
            ellist.append(bb)
            filtlist.append(cc)
    
    if freq == 148:
        ellist = []
        filtlist = []
        files = open('/u/sudeep/ACT/mapAnalysis/compactSources/data/B_l_AR2.dat','r')
        for line in files:
            aa = line.split(' ')
            bb = float(aa[0])
            cc = 1./float(aa[1])
            ellist.append(bb)
            filtlist.append(cc)

    kMap = ftthing.kMap.copy()
    kFilter0 = numpy.real(kMap.copy())*0.+ 1.
    kFilter = kMap.copy()*0.
    FlSpline = splrep(ellist,filtlist,k=3)
    ll = numpy.ravel(ft148.modLMap)
    kk = (splev(ll,FlSpline))
    kFilter = numpy.reshape(kk,[ftthing.Ny,ftthing.Nx])
    kFilter0 *= kFilter
    kMap[:,:] = kFilter0[:,:]

    return kMap

def deconvolveMap(conmap,fwhmRad):

    #deconvolving####################we are only using ls up to 20000 to prevent leakage
    ell = numpy.arange(20000)
    Filter = numpy.exp(ell**2.0*fwhmRad**2.0/(16.*numpy.log(2.)))
    print "deconvolving.."
    ft2 = fftTools.fftFromLiteMap(conmap)
    kMap = ft2.kMap.copy()
    kFilter0 = numpy.real(kMap.copy())*0.+ 1.
    kFilter = kMap.copy()*0.
    FlSpline = splrep(ell,Filter,k=3)
    ll = numpy.ravel(ft2.modLMap)
    kk = (splev(ll,FlSpline))
    kFilter = numpy.reshape(kk,[conmap.Ny,conmap.Nx])
    kFilter0 *= kFilter
    kMap[:,:] *= kFilter0[:,:]
    ft2.kMap = kMap
    deconlensedMap2data = numpy.real(ifft2(kMap))
    conmap.data[:] = deconlensedMap2data[:]
    return conmap


def createKMap(map0, map1, LF):
    
    LF.map1 = map0.copy()
    LF.map2 = map1.copy()

    lpMap  = LF.applyLowPass()
    hpMap  = LF.applyHighPass()

    ftkap, kappaMap, norma = LF.makeKappaMap(trimAtL=trimNofLAtL)
    
    LF.map1 = map1.copy()
    LF.map2 = map0.copy()
    
    lpMap  = LF.applyLowPass()
    hpMap  = LF.applyHighPass()

    ftkap, kappaMap2, norma = LF.makeKappaMap(trimAtL=trimNofLAtL)
    
    #bla = ftkap.copy()
    #bla.powerMap = numpy.fft.ifftshift(norma)
    #lL,lU,lBin,clN,err,w = bla.binInAnnuli(os.environ['LENSRECONS_DIR']+os.path.sep+'params/BIN_100_LOG')

    return ftkap, (kappaMap+kappaMap2)*LF.powerMap1.pixScaleX*LF.powerMap1.pixScaleY/2.#map1.pixScaleX*map1.pixScaleY/2.

def createKMapOne(map0, LF):
    
    LF.map1 = map0.copy()
    LF.map2 = map0.copy()

    lpMap  = LF.applyLowPass()
    hpMap  = LF.applyHighPass()

    ftkap, kappaMap, norma, noiseBiasNew = LF.makeKappaMap(trimAtL=trimNofLAtL)
    
    #bla = ftkap.copy()
    #bla.powerMap = numpy.fft.ifftshift(norma)
    #lL,lU,lBin,clN,err,w = bla.binInAnnuli(os.environ['LENSRECONS_DIR']+os.path.sep+'params/BIN_100_LOG')

    return ftkap, (kappaMap)*LF.powerMap1.pixScaleX*LF.powerMap1.pixScaleY, norma, noiseBiasNew

def createDivMap(map0, map1, LF):
    '''
    @brief Given two maps creates the divergence of TGradT map. The high and low pass
    filters are symmetrized between the maps.
    @param map0  the first map
    @param map1 the second map
    @param LF lensFilter instanace 
    '''

    LF.map0 = map0.copy()
    LF.map1 = map1.copy()

    lpMap  = LF.applyLowPass()
    hpMap  = LF.applyHighPass()

    gradLP = lpMap.takeGradient()

    hpGradLpX = lpMap.copy()
    hpGradLpY = lpMap.copy()

    hpGradLpX.data[:] = gradLP.gradX.data[:]*hpMap.data[:]
    hpGradLpY.data[:] = gradLP.gradY.data[:]*hpMap.data[:]

    divMap0 = liteMap.takeDivergence(hpGradLpX,hpGradLpY)

    LF.map0 = map1.copy()
    LF.map1 = map0.copy()

    lpMap  = LF.applyLowPass()
    hpMap  = LF.applyHighPass()

    gradLP = lpMap.takeGradient()

    hpGradLpX = lpMap.copy()
    hpGradLpY = lpMap.copy()

    hpGradLpX.data[:] = gradLP.gradX.data[:]*hpMap.data[:]
    hpGradLpY.data[:] = gradLP.gradY.data[:]*hpMap.data[:]

    divMap1 = liteMap.takeDivergence(hpGradLpX,hpGradLpY)

    divMapMean = divMap1.copy()

    divMapMean.data[:] = (divMap1.data[:]+ divMap0.data[:])/2.

#BLAKE: get the divmap1 by flipping map0 and map1, and return mean

    return divMapMean #should be divMapMean

def extrapolateToPowerMap(crossPowerN,ell,Filter):

     #ft2 = fftTools.fftFromLiteMap(conmap)
    kMap = crossPowerN.powerMap.copy()#ft2.kMap.copy()
    kFilter0 = numpy.real(kMap.copy())*0.+ 1.
    kFilter = kMap.copy()*0.
    FlSpline = splrep(ell,Filter,k=3)
    ll = numpy.ravel(crossPowerN.modLMap)
    kk = (splev(ll,FlSpline))
    kFilter = numpy.reshape(kk,[crossPowerN.Ny,crossPowerN.Nx])
    kFilter0 *= kFilter
    kMap[:,:] *= kFilter0[:,:]*0.
    kMap[:,:] += kFilter0[:,:]
    copy = crossPowerN.copy()
    copy.powerMap = kMap
    return copy
########################################################

p = flipperDict.flipperDict()
p.read_from_file(sys.argv[1])
resultDir = p['resultDir']
nMC = p['nMC']
os.system("cp /u/bsherwin/codes/lensRecons/bin/mcnew %s/mcnewback"%resultDir)
os.system("cp /u/bsherwin/codes/lensRecons/python/lensT.py %s/lensTBackup"%resultDir)

for q in xrange(1):#nMC):
    print 'q=', q
    tMapFile = p['tMapFile']#'/scr/queequeg1/shared/projects/lensing/bcsCrossCorr/act148Shifted/patchesOnBCSFilteredShifted0/patch000_4'#+ str(q)#dir where maps are stored#A #list of mapNames
    tMapFile2 = p['tMapFile220']
    kappaMapFile = p['kappaFile']#'/scr/queequeg1/shared/projects/lensing/mcSims/seasonNoise1.0/kappa.fits'
    gMapFile = p['gMapFile']#'/scr/queequeg1/shared/projects/lensing/bcsCrossCorr/bcsMaps/kappaMap/patch000_0' #'/scr/queequeg1/shared/projects/#lensing/mcSims/seasonNoise17.0/patches_lensed2D/#patchx'
    accuracyBoost = p['accuracyBoost']
    apodize = p['apodize']
    resultDir = p['resultDir']
    noiseUkArcmin = p['noiseUkArcmin']
    thetaFWHMarcmin = p['beamFWHM']
    cosSqCutoffParams=p['filterCosSqCutoff']
    trimNofLAtL = p['trimNofLAtL']
    lMaxForPlots = p['lMaxForPlots']
    plotDiagnostics = p['plotDiagnostics']
    TCMB = 2.726e6
    try:
        os.mkdir(resultDir)
    except:
        pass

    pw = prewhitener.prewhitener(1., 0.02, 0.)
    

    thetaFWHM = thetaFWHMarcmin*numpy.pi/(180.*60.)

    #True kappa Map
    kappaTrue = liteMap.liteMapFromFits(kappaMapFile)
    
    
    #true auto spectrum
    p2d = fftTools.powerFromLiteMap(kappaTrue,applySlepianTaper=True)
    lL,lU,lBin,clKappaTrueBin,err,w = p2d.binInAnnuli('/u/bsherwin/flipper/params/BIN_100_LOG')

    #read in the theory spectra: unlensed Cls and clKappaKappa
    print "Reading in theory spectra"
    X = pylab.load(os.environ['LENSRECONS_DIR']+os.path.sep+'params/wmap5scalCls.dat')
    l = X[:,0]
    lsqClOver2pi = X[:,1]
    cl = lsqClOver2pi*(2*numpy.pi)/(l*(l+1.0))
    cl /=TCMB**2 #unlensed cls
    clkap = X[:,4]/(4.*TCMB**2) #clkapkap

    X = pylab.load(os.environ['LENSRECONS_DIR']+os.path.sep+'params/wmap5lensedCls.dat')
    lLensed = X[:,0]
    lsqClOver2pi = X[:,1]
    clLensed = lsqClOver2pi*(2*numpy.pi)/(lLensed*(lLensed+1.0))
    clLensed /=TCMB**2 #lensed cls

    transferFl = pw.correctSpectrum(l,(l*0.+1.))
    transferFlL = pw.correctSpectrum(lLensed,(lLensed*0.+1.))
    #cl = cl/transferFl
    #clLensed = clLensed/transferFlL
    
    if noiseUkArcmin == None:
        deltaT = 0.0
    else:
        deltaT = noiseUkArcmin/thetaFWHMarcmin

    #the theoretical noise spectrum
    nl = (deltaT*thetaFWHM)**2*numpy.exp(l*(l+1)*thetaFWHM**2/(8.*numpy.log(2.)))/TCMB**2
    
    ellist = []
    gktheory = []
    ckktheory = []
    nkktheory = []
    deltatheory = []
    files = open('/scr/queequeg1/shared/projects/lensing/bcsCrossCorr/bcsMaps/analyticCalc/clgk_BCS0.1to1.5_.dat','r')
    for line in files:
        aa = line.split(' ')
        bb = float(aa[0])
        cc = float(aa[1])
        ee = float(aa[3])
        dd = float(aa[4])
        ff = float(aa[2])/10.
        ellist.append(bb)
        gktheory.append(cc)
        ckktheory.append(ee)
        nkktheory.append(dd)
        deltatheory.append(ff)
    ellist2 = []
    ggtheory = []
    files = open('/scr/queequeg1/shared/projects/lensing/bcsCrossCorr/bcsMaps/analyticCalc/clgg_BCS0.1to1.5_.dat','r')
    for line in files:
        aa = line.split(' ')
        bb = float(aa[0])
        cc = float(aa[1])
        dd = float(aa[2])
        ellist2.append(bb)
        ggtheory.append(cc+dd)
    
    bla = nkktheory

    for i in xrange(len(bla)):
        bla[i] += ckktheory[i]
    
    print "lengths", len(ellist),len(ellist2),len(gktheory),len(ggtheory),len(ckktheory),len(nkktheory),len(bla)
    

    #nMaps = len(mapNames)

    map = liteMap.liteMapFromFits(tMapFile)#+'patch000_4')
    mask = map.copy()
    del map
    if apodize:
        print "Will apodize edges..."
        mask.data[:] = 0.
        mask.data[10:mask.Ny-10,10:mask.Nx-10] = 1.
        mask = mask.convolveWithGaussian(fwhm=10.)
        #taper = utils.slepianTaper00(mask.Nx, mask.Ny, 2.0)
    else:
        mask.data[:] = 1.0
    value = mask.data.copy()
    value[:] = 1.0
    sumorig = value.sum()
    sumnew = mask.data.sum()
    fraction = sumnew/sumorig
    print 'map apod fraction', fraction
    taper = utils.slepianTaper00(mask.Nx,mask.Ny,3.0)
    taper[:] = 1.
    mask.data[:] /= numpy.sqrt(fraction)  
    
            
    #kappaTrue.data[:] *= mask.data[:]
    p2d = fftTools.powerFromLiteMap(kappaTrue,applySlepianTaper=False)
    p22d = p2d.trimAtL(trimNofLAtL)
    lL,lU,lBin,clKappaTrueBin,err,w = p22d.binInAnnuli('/u/bsherwin/flipper/params/BIN_100_LOG')

    irfactor = 1.#0.*10.*10.
    #BLAKE read in other necessary parameters, theory etc
    mapp = liteMap.liteMapFromFits(tMapFile)#+'patch000_4')
    #mapp.convertToMicroKFromJyPerSr(148., Tcmb = 2.726)

    mapt1 = liteMap.liteMapFromFits(tMapFile)#+'patch000_0')
    #mapt1.convertToMicroKFromJyPerSr(148., Tcmb = 2.726)
    mapt2 = liteMap.liteMapFromFits(tMapFile)#+'patch000_1')
    #mapt2.convertToMicroKFromJyPerSr(148., Tcmb = 2.726)
    mapt3 = liteMap.liteMapFromFits(tMapFile)#+'patch000_2')
    #mapt3.convertToMicroKFromJyPerSr(148., Tcmb = 2.726)
    mapt4 = liteMap.liteMapFromFits(tMapFile)#+'patch000_3')
    #mapt4.convertToMicroKFromJyPerSr(148., Tcmb = 2.726)
    mapt1.data[:] /= TCMB/irfactor#*0.83#*taper[:]
    mapt2.data[:] /= TCMB/irfactor#*0.83#*taper[:]
    mapt3.data[:] /= TCMB/irfactor#*0.83#*taper[:]
    mapt4.data[:] /= TCMB/irfactor#*0.83#*taper[:]

    #mapnull = mapp.copy()
    #mapnull.data[:] = (mapt1.data[:]-mapt2.data[:])/2./numpy.sqrt(2.)
    mapsplit1 = mapp.copy()
    mapsplit2 = mapp.copy()
    mapsplit1.data[:] = (mapt1.data[:] + mapt2.data[:])/2.
    mapsplit2.data[:] = (mapt3.data[:] + mapt4.data[:])/2.

    
    a220mapp = liteMap.liteMapFromFits(tMapFile2)#+'patch000_4')
    a220mapt1 = liteMap.liteMapFromFits(tMapFile2)#+'patch000_0')
    a220mapt2 = liteMap.liteMapFromFits(tMapFile2)#+'patch000_1')
    a220mapt3 = liteMap.liteMapFromFits(tMapFile2)#+'patch000_2')
    a220mapt4 = liteMap.liteMapFromFits(tMapFile2)#+'patch000_3')
    a220mapt1.data[:] *= mask.data[:]/TCMB#*taper[:]
    a220mapt2.data[:] *= mask.data[:]/TCMB#*taper[:]
    a220mapt3.data[:] *= mask.data[:]/TCMB#*taper[:]
    a220mapt4.data[:] *= mask.data[:]/TCMB#*taper[:]
    a220mapp.data[:] -= a220mapp.data.mean()
    a220mapt1.data[:] -= a220mapt1.data.mean()
    a220mapt2.data[:] -= a220mapt2.data.mean()
    a220mapt3.data[:] -= a220mapt3.data.mean()
    a220mapt4.data[:] -= a220mapt4.data.mean()
    a220mapnull = a220mapp.copy()
    a220mapnull.data[:] = (a220mapt1.data[:]-a220mapt2.data[:])/2./numpy.sqrt(2.)
    a220mapsplit1 = a220mapp.copy()
    a220mapsplit2 = a220mapp.copy()
    a220mapsplit1.data[:] = (a220mapt1.data[:] + a220mapt2.data[:])/2.
    a220mapsplit2.data[:] = (a220mapt3.data[:] + a220mapt4.data[:])/2.
    
        #mapp2 = deconvolveMap(mapp,thetaFWHM)
    mapp.data[:] /= TCMB/irfactor#*0.83#/0.8524**4.
    a220mapp.data[:] = (a220mapsplit1.data[:] + a220mapsplit2.data[:])/2.#*= mask.data[:]/TCMB#/0.8524**4.

    thetaFWHM = thetaFWHMarcmin*numpy.pi/(180.*60.)
    deltaT = noiseUkArcmin/thetaFWHMarcmin
    pixSide = numpy.sqrt(mapp.pixScaleX*mapp.pixScaleY)*180.*60./numpy.pi
    deltaTPerPix = noiseUkArcmin/(pixSide) #noise per sky pixel

#generate and add Gaussian noise

    #mapp.data[:] = (deltaTPerPix*numpy.random.randn(mapp.Ny,mapp.Nx)/TCMB)[:]
    #mapt1.data[:] = (deltaTPerPix*numpy.random.randn(mapp.Ny,mapp.Nx)/TCMB)[:]
    #mapt2.data[:] = (deltaTPerPix*numpy.random.randn(mapp.Ny,mapp.Nx)/TCMB)[:]
    #mapt3.data[:] = (deltaTPerPix*numpy.random.randn(mapp.Ny,mapp.Nx)/TCMB)[:]
    #mapt4.data[:] = (deltaTPerPix*numpy.random.randn(mapp.Ny,mapp.Nx)/TCMB)[:]
    mapp.data[:] *= mask.data[:]
    mapt1.data[:] *= mask.data[:]
    mapt2.data[:] *= mask.data[:]
    mapt3.data[:] *= mask.data[:]
    mapt4.data[:] *= mask.data[:]
    mapp.data[:] -= mapp.data.mean()
    mapt1.data[:] -= mapt1.data.mean()
    mapt2.data[:] -= mapt2.data.mean()
    mapt3.data[:] -= mapt3.data.mean()
    mapt4.data[:] -= mapt4.data.mean()

    mapsAll = [mapt1,mapt2,mapt3,mapt4]
    a220mapsAll = [a220mapt1,a220mapt2,a220mapt3,a220mapt4]
    themapppower = fftTools.powerFromLiteMap(mapp)
    a220themapppower = fftTools.powerFromLiteMap(a220mapp)
    
    lBin,clBinCrossMean,meanCrossPM = meanCrossSpec(mapsAll,applySlepianTaper=False)
    lBin,clBinAutoMean,meanAutoPM = meanAutoSpec(mapsAll,applySlepianTaper=False)

    pylab.clf()
    pylab.loglog(lBin,clBinAutoMean)
    pylab.loglog(l,cl)
    pylab.savefig('checkIRMap.png')
    pylab.clf()
    
    #noiseTh2D = (themapppower.powerMap - meanCrossPM)#/2./numpy.sqrt(2.)(meanAutoPM-meanCrossPM)/4.#
    #noiseTh2D = (meanAutoPM-meanCrossPM)/4.
    #noiseThDiff = - noiseTh2D2 + noiseThAlt
    #noiseTh2D = noiseThDiff
    #get mean cross spectrum - compare with theory
    irCross = extrapolateToPowerMap(themapppower,lBin,clBinCrossMean)
    irNoise = extrapolateToPowerMap(themapppower,l,nl)

    lBin,a220clBinCrossMean,a220meanCrossPM = meanCrossSpec(a220mapsAll,applySlepianTaper=False)
    lBin,a220clBinCrossMean,a220meanAutoPM = meanAutoSpec(a220mapsAll,applySlepianTaper=False)
    a220noiseTh2D = (a220meanAutoPM - a220meanCrossPM)/4.####(a220meanAutoPM-a220meanCrossPM)/4.#########/2.#######/numpy.sqrt(2.)

    ft148 = fftTools.fftFromLiteMap(mapp)
    ft220 = fftTools.fftFromLiteMap(a220mapp)
    new148power = fftTools.powerFromLiteMap(mapp)
    
    
    clkMap = ft148.kMap.copy()
    kFilter0 = numpy.real(clkMap.copy())*0.+ 1.
    kFilter = clkMap.copy()*0.
    FlSpline = splrep(l,cl,k=3)
    ll = numpy.ravel(ft148.modLMap)
    kk = (splev(ll,FlSpline))
    kFilter = numpy.reshape(kk,[ft148.Ny,ft148.Nx])
    kFilter0 *= kFilter
    clkMap[:,:] = kFilter0[:,:]

    fttot = ft148.copy()
    #fttot.kMap[:,:] *= factor148[:,:]
    fwhmRad = thetaFWHMarcmin*numpy.pi/(180.*60.)
    ell = numpy.arange(20000)
    FilterBeam = numpy.exp(ell**2.0*fwhmRad**2.0/(16.*numpy.log(2.)))
    beamPower = extrapolateToPowerMap(themapppower,ell,FilterBeam)
    a148Beam = beamPower.powerMap#deconvolveMapRealBeam(148,fttot)
    #a220Beam = 1.+0.*noiseTh2D[:,:]#deconvolveMapRealBeam(148,fttot)
    #factor148 = fttot.kMap.copy()
    #factor220 = fttot.kMap.copy()
    #factor148[:,:] = 1./(clkMap[:,:] + noiseTh2D[:,:]*a148Beam[:,:]**2.)
    #factor220[:,:] = 1./(clkMap[:,:] + a220noiseTh2D[:,:]*a220Beam[:,:]**2.)
    #fttot.kMap[:,:] = ft148.kMap[:,:]*a148Beam[:,:]*factor148[:,:] + ft220.kMap[:,:]*a220Beam[:,:]*factor220[:,:]
    #fttot.kMap[:,:] *= 1./(factor148[:,:]+factor220[:,:])
    #fttot.kMap[:,:] *= a148Beam[:,:]
    powerOfComb = powerFromKMap(fttot.kMap,fttot.kMap,themapppower,themapppower)
    
    #noiseTotal = noiseTh2D.copy()
    #noiseTotal[:,:] = noiseTh2D[:,:]*a148Beam[:,:]**2.*(factor148[:,:])**2. + a220noiseTh2D*a220Beam[:,:]**2.*(factor220[:,:])**2.
    #noiseTotal[:,:] *= 1./(factor148[:,:]+factor220[:,:])**2.
    #noiseTotal[:,:] = noiseTh2D[:,:]*a148Beam[:,:]**2.####TEMPORARY!!!

    pylab.clf()
    holder = p2d.copy()
    #holder.powerMap = noiseTotal
    #holder.plot(zoomUptoL=3000,log=True)
    #pylab.savefig('2DnoiseDiff.png')
    #pylab.clf()
    
    #sumMap = holder.copy()
    #sumMap.powerMap[:,:] = holder.powerMap[:,:] + clkMap[:,:]
    #new148power.powerMap[:,:] *= a148Beam[:,:]**2.
    #crossMap = holder.copy()
    #crossMap.powerMap[:,:] = meanCrossPM[:,:]



 
    #lL,lU,lBin,a148P, plErr_12,weights = new148power.binInAnnuli('/u/bsherwin/flipper/params/BIN_100_LOG')
    #lL,lU,lBin,sumP, plErr_12,weights = sumMap.binInAnnuli('/u/bsherwin/flipper/params/BIN_100_LOG')
    lL,lU,lBin,noiseP, plErr_12,weights = holder.binInAnnuli('/u/bsherwin/flipper/params/BIN_100_LOG')
    lL,lU,lBin,mapP, plErr_12,weights = powerOfComb.binInAnnuli('/u/bsherwin/flipper/params/BIN_100_LOG')
    #lL,lU,lBin,crossP, plErr_12,weights = crossMap.binInAnnuli('/u/bsherwin/flipper/params/BIN_100_LOG')
    transferF = pw.correctSpectrum(lBin,(lBin*0.+1.))


    #pylab.loglog(lBin,transferF*a148P*lBin**2.,label='148P')
    pylab.loglog(lBin,noiseP*lBin**2.,label='noiseP')
    pylab.loglog(lBin,mapP*lBin**2.,label='mapP')
    #pylab.loglog(lBin,transferF*sumP*lBin**2.,label='sumP')
    #pylab.loglog(lBin,transferF*crossP*lBin**2.,label='crossP')
    pylab.loglog(l,cl*l**2.,label='cl')
    pylab.legend(loc='best')
    pylab.savefig('noisePlots2.png')
    pylab.clf()

    print "now"

    transferFlN = pw.correctSpectrum(l,(l*0.+1.))
    trPower = extrapolateToPowerMap(holder,l,transferFlN)
    trPowerSq = extrapolateToPowerMap(holder,l,numpy.sqrt(transferFlN))
    tr = trPower.powerMap.copy()
    trSq = trPowerSq.powerMap.copy()
    
    #powerOfComb.powerMap[:,:] *= tr[:,:]

    
    #get mean cross spectrum - compare with theory
    #ft148.kMap[:,:]*a148Beam[:,:]
    #initialize LF (for now pass in the theory for all spectra)
    LF = lensFiltersElSpace(fttot.kMap\
                            ,fttot.kMap\
                            ,powerOfComb\
                            ,l\
                            ,cl\
                            ,irNoise.powerMap*0.\
                            ,lBin\
                            #,clBinCrossMean\
                            #,clBinCrossMean\
                            #,clBinCrossMean\
                            #,thetaFWHM\
                            #,cosSqCutoffParams)
                            ,irCross.powerMap\
                            ,irCross.powerMap\
                            ,irCross.powerMap\
                            ,irCross.powerMap\
                            ,thetaFWHM\
                            ,cosSqCutoffParams)#themappower replaces meanCrossPM
                    
 

    
    divMapsAll = []
    indicesList = []

    #LF.map1 = mapsplit1
    #LF.map2 = mapsplit2
    #bla = LF.applyHighPass()
    #bla = LF.applyLowPass()
   # LF.normalizeTwoD(trimAtL = trimNofLAtL)
   # pickle.dump(LF.normalization[0:6000],open('normalizationaltorig.pkl','w'))
   # pickle.dump(LF.noiseBias[0:6000],open('noiseBiasaltorig.pkl','w'))
    #sys.exit()

   
    ftex, ftkappaMap, normOutput, noiseBiasNew = createKMapOne(fttot.kMap, LF)
    #ftex, ftkappaMap = createKMap(fttot.kMap*trSq,fttot.kMap*trSq,LF)
    #ftex2, ftkappaMapNull = createKMap(mapsplit1,mapsplit2, LF2) 
    #ftex2, ftkappaMapNull = createKMap(mapnull,mapnull,LF)
    gal = liteMap.liteMapFromFits(kappaMapFile)
    gal.data[:] *= mask.data[:]
    galf = fftTools.fftFromLiteMap(gal)
    galMap = galf.kMap.copy()
    galMapT = trimShiftKMap(galMap,trimNofLAtL,0,0,galf.lx,galf.ly)
    #galpower0 = powerFromKMap(galMapT,galMapT,p2d,ftex)
    #lL,lU,lBin,divfactor, plErr_12,weights = galpower0.binInAnnuli('/u/bsherwin/flipper/params/BIN_100_LOG')  

    #kMap = ftex.powerMap.copy()
    #kFilter0 = numpy.real(kMap.copy())*0.+ 1.
    #kFilter = kMap.copy()*0.
    #FlSpline = splrep(lBin,divfactor,k=3)
    #ll = numpy.ravel(ftex.modLMap)
    #kk = (splev(ll,FlSpline))
    #kFilter = numpy.reshape(kk,[ftex.Ny,ftex.Nx])
    #kFilter0 *= kFilter
    #kMap[:,:] = kFilter0[:,:]
    bla = ftex.copy()
    bla.powerMap = noiseBiasNew


    #normOutput = noiseBiasNew
    #nB = powerFromKMap(normOutput,normOutput,p2d,ftex)#ftex.copy()
    #nB.powerMap = noiseBiasNew
    #nB.powerMap = normOutput
    
    #galMapT[:,:] /= kMap[:,:]#galpower0.powerMap[:,:]
    crossPower0 = powerFromKMap(ftkappaMap,galMapT,p2d,ftex)
    #crossPowerNull = powerFromKMap(ftkappaMapNull,galMapT,p2d,ftex2)
    scatter = powerFromKMap(ftkappaMap,ftkappaMap,p2d,ftex)
   # galpower = powerFromKMap(galMapT,galMapT,p2d,ftex)
    #nullft = powerFromKMap(ftkappaMapNull,ftkappaMapNull,p2d,ftex2)
    #error = scatter.copy()
    #error.powerMap[:] = numpy.sqrt(numpy.abs(scatter.powerMap[:]*galpower.powerMap[:]))
    
    #pylab.clf()
    #crossPower0.plot()
    #pylab.savefig('nullPower.png')
    #pylab.clf()

    #scatter.powerMap[:] *= galpower.powerMap[:]
    #scatter.powerMap[:] = numpy.sqrt(numpy.abs(scatter.powerMap[:]))
    
    #ftex, ftKappaMap1 = createKMap(mapsAll[0], mapsAll[1], LF)
    #ftex, ftKappaMap2 = createKMap(mapsAll[2], mapsAll[3], LF)
    #divMapsAll += [ftKappaMap1]
    #divMapsAll += [ftKappaMap2]
    #indicesList += [[0,1]]
    #indicesList += [[2,3]]
    #biases = []




  #  crossPower0 = fftTools.powerFromLiteMap(divMapsAll[0],divMapsAll[0],applySlepianTaper=False)
    lL,lU,lBin,kgal, plErr_12,weights = crossPower0.binInAnnuli('/u/bsherwin/flipper/params/BIN_100_LOG')
    #lL,lU,lBin,knullgal, plErr_12,weights = crossPowerNull.binInAnnuli(os.environ['LENSRECONS_DIR']+os.path.sep+'params/BIN_100_LOG')
    lL,lU,lBin,kk, plErr_12,weights = scatter.binInAnnuli('/u/bsherwin/flipper/params/BIN_100_LOG')
    lL,lU,lBin,newNoiseBias, plErr_12,weights = bla.binInAnnuli('/u/bsherwin/flipper/params/BIN_100_LOG')
    #lL,lU,lBin,galgal, plErr_12,weights = galpower.binInAnnuli('/u/bsherwin/flipper/params/BIN_100_LOG')
    #lL,lU,lBin,knullknull, plErr_12,weights = nullft.binInAnnuli(os.environ['LENSRECONS_DIR']+os.path.sep+'params/BIN_100_LOG')
    #lL,lU,lBin,err, plErr_12,weights = error.binInAnnuli('/u/bsherwin/flipper/params/BIN_100_LOG')
    #lL,lU,lBin,err, plErr_12,weights = error.binInAnnuli('/u/bsherwin/flipper/params/BIN_100_LOG')
    #lL,lU,lBin,nBB, plErr_12,weights = nB.binInAnnuli('/u/bsherwin/flipper/params/BIN_100_LOG')

    #err /= numpy.sqrt(weights/2.)

    #print "product is:", product
    #print "scatter is:", scatter
    #sys.exit()
    pylab.clf()
   
    #pylab.loglog(lBin,clKappaTrueBin,label='true map auto')
    thefactorl = numpy.sqrt((2.*numpy.pi)/2./lBin/100./mapp.Nx/mapp.Ny/mapp.pixScaleX/mapp.pixScaleY)

    pylab.loglog(lBin,kk,label='reconstructed kappa')
    pylab.loglog(lBin,-kgal,'o',label='cross')
    #pylab.loglog(lBin,kgal,'o',label='cross')
    #pylab.loglog(lBin,numpy.abs(productnull)+sigma,label='null + Delta')
    #pylab.loglog(lBin,numpy.abs(productnull)-sigma,label='null - Delta')
    #pylab.loglog(lBin,clNoiseBias,label='noise bias')
    pylab.loglog(lBin,clKappaTrueBin,label='theory kappa')
    #pylab.loglog(lBin,clKappaTrueBin+normalize,label='true + normalization')
    #pylab.loglog(lBin,product-clNoiseBias,label='crosstmap - noise bias')
    #pylab.loglog(lBin,clKappaTrueBin+noiseBias,label='true + noise bias -nonBin')
    #pylab.loglog(lBin,product-noiseBias,label='crosstmap - noise bias-nonBin')
 
    pylab.xlabel(r'$\ell$')
    pylab.ylabel(r'$C_\ell$')
    pylab.xlim(10,lMaxForPlots)
    #pylab.ylim(10.**(-10.),10.**(-4.))
    pylab.legend(loc='best')
    #utils.saveAndShow()
    pylab.savefig(resultDir+os.path.sep+'powerSpectraKappa.png')
    pylab.clf()

    pickle.dump(lBin,open('lBin2.pkl','w'))
    pickle.dump(thefactorl,open('factor2.pkl','w'))
    pickle.dump(kk,open('kk2.pkl','w'))
    pickle.dump(clKappaTrueBin,open('clKappaTrue2.pkl','w'))
    pickle.dump(newNoiseBias,open('newNoiseBias.pkl','w'))
    #pickle.dump(noiseBiasNew,open('origNB.pkl','w'))
#pickle.dump(err,open('err.pkl','w'))

    #pylab.plot(lBin,-knullgal,label='knullgal')
   # pylab.plot(lBin,-kgal,label='kgal')
    #pylab.plot(lBin,err,label='err')
    #pylab.plot(lBin,-err,label='err')
    #pylab.plot(lBin,-kgal,label='kgal')
    #pylab.xlim(10,lMaxForPlots)
    #pylab.legend(loc='best')
    #utils.saveAndShow()
   # pylab.savefig(resultDir+os.path.sep+'2ndpsk.png')
    #pylab.clf()

    #pylab.loglog(lBin,clKappaTrueBin,label='clkappa')
    #pylab.loglog(lBin,Delta,label='DeltaCl')
    #pylab.legend(loc='best')
    #pylab.savefig(resultDir+os.path.sep+'signalnoise.png')
    #pylab.clf()

    #pylab.loglog(lLensed,clkap,label='clkap')
    #pylab.loglog(lLensed,lLensed*lLensed*nl/840.,label='noise/840')
    #pylab.legend(loc='best')
    #pylab.savefig(resultDir+os.path.sep+'doweneed8pt.png')
    #pylab.clf()

   


