Commit 9cef3e5c authored by Maude Le Jeune's avatar Maude Le Jeune
Browse files

ps subtraction

parent 061d2da4
......@@ -31,8 +31,7 @@
using namespace std;
void gaussfit (Healpix_Map<double> &hmap, Healpix_Map<double> &smap, double* tabcat, int ns,int ncol, double fwhm_arcmin,double nfwhm ,double latcut,double maxxi2 ) {
void gaussfit (Healpix_Map<double> &hmap, Healpix_Map<double> &smap, double* tabcat, int ns,int ncol, double fwhm_arcmin,double nfwhm ,double latcut,double maxxi2 , int nofit) {
bool sigmap = true;
......@@ -50,9 +49,9 @@ void gaussfit (Healpix_Map<double> &hmap, Healpix_Map<double> &smap, double* tab
struct timeval startTime;
struct timeval endTime;
#pragma omp parallel
{
#pragma omp for
//#pragma omp parallel
// {
//#pragma omp for
for (int i=0; i<ns; i++){
// for(s = cat.begin(); s != cat.end(); s++){
// pointing p = (*s).getPointing();
......@@ -60,18 +59,34 @@ void gaussfit (Healpix_Map<double> &hmap, Healpix_Map<double> &smap, double* tab
double guess[3];
double err[3];
pointing p;
p.phi = tabcat[icol+1]* degr2rad;
p.theta = tabcat[icol+2];
p.theta = (90. - p.theta) * degr2rad;
guess[1] = p.theta;
guess[2] = p.phi;
guess[0] = tabcat[icol+3];
if (abs(90-(p.theta*180.0/pi)) < latcut)
chisqrdof = INFINITY;
else {
healpixGaussianFit(hmap, fwhm, guess, err, &chisqrdof,smap,radius);
double ei;
if (nofit==0){
p.phi = tabcat[icol+1]* degr2rad;
p.theta = tabcat[icol+2];
p.theta = (90. - p.theta) * degr2rad;
guess[1] = p.theta;
guess[2] = p.phi;
guess[0] = tabcat[icol+3];
ei = tabcat[icol+4];
if (abs(90-(p.theta*180.0/pi)) < latcut)
chisqrdof = INFINITY;
else {
healpixGaussianFit(hmap, fwhm, guess, err, &chisqrdof,smap,radius);
}
}
else {
guess[1] = tabcat[icol+0] ;
guess[2] = tabcat[icol+1];
guess[0] = tabcat[icol+2] ;
err[1] = tabcat[icol+3];
err[2] = tabcat[icol+4] ;
err[0] = tabcat[icol+5] ;
chisqrdof = tabcat[icol+6];
}
// fit is ok, set to new value
if (chisqrdof < maxxi2 && guess[0]>0 && err[0]!=INFINITY && err[1]<guess[1]){
addGaussianPatternOnSky(hmap, pointing(guess[1],guess[2]), fwhm, guess[0], err[0]/10, true);
......@@ -85,16 +100,18 @@ void gaussfit (Healpix_Map<double> &hmap, Healpix_Map<double> &smap, double* tab
}
// fit is bad, keep old values
else {
else if (nofit==0) {
tabcat[icol+0] = p.theta;
tabcat[icol+1] = p.phi;
tabcat[icol+2] = tabcat[icol+3];
tabcat[icol+3] = 0;
tabcat[icol+5] = tabcat[icol+4];
tabcat[icol+4] = 0;
if (ei==0)
printf("Warning: null value\n");
tabcat[icol+5] = ei;
tabcat[icol+6] = chisqrdof;
}
//i++;
}
// }
}
}
......@@ -192,7 +192,7 @@ def write_table (filename, mat, colname=None):
if colname is None:
colname = []
for c in range(ncols):
colname.append("column "+str(c))
colname.append("c"+str(c))
cols = []
for c in range(ncols):
......
......@@ -178,9 +178,30 @@ def ps_gaussfit (map, smap, cat, fwhm_arcmin, nsigma=1, latcut=0, maxxi2=10):
nsource = shape(cat)[0]
Bcat = zeros((nsource, 7))
Bcat[:,0:5] = cat
spherelib_map._psgaussfit(map, smap, Bcat, fwhm_arcmin, nsigma, latcut, maxxi2)
spherelib_map._psgaussfit(map, smap, Bcat, fwhm_arcmin, nsigma, latcut, maxxi2, 0)
return Bcat
def ps_subtract (map, cat, fwhm_arcmin, nsigma=1, latcut=0, maxxi2=10):
""" Perform gaussian fit.
Parameters
----------
map: healpix map
smap: healpix map of the background
cat : ps catalog array like (nsource, 5) with columns corresponding to (pixel, glon, glat, Ampl, sigmaAmpl)
fwhm_arcmin: map resolution
nsigma: fit radius
latcut: cut on latitude (degree)
maxxi2: maximum x2 value
Returns
-------
array-like (nsource, 7) with columns corresponding to (theta, phi, Ampl, s_theta, s_phi, s_flux, xi2/dof)
"""
smap = zeros(shape(map))
spherelib_map._psgaussfit(map, smap, cat, fwhm_arcmin, nsigma, latcut, maxxi2, 1)
def cat2mask (mask, theta, phi, fwhm, snr=None):
""" Compute an Healpix map with null values around the
sources.
......
......@@ -357,7 +357,7 @@ void _psdetection(double* tabmap, int npix, double** tabcat, int* ncat, double f
*ncat = ns*ncol;
}
void _psgaussfit (double* tabmap, int npix, double* sigmamap, int npix2, double* tabcat, int ns,int ncol, double fwhm_arcmin,double nfwhm ,double latcut,double maxxi2 ) {
void _psgaussfit (double* tabmap, int npix, double* sigmamap, int npix2, double* tabcat, int ns,int ncol, double fwhm_arcmin,double nfwhm ,double latcut,double maxxi2, int nofit ) {
if (npix != npix2)
PyErr_Format(PyExc_ValueError,"Arrays of different lengths (%d %d)",npix, npix2);
......@@ -370,7 +370,7 @@ void _psgaussfit (double* tabmap, int npix, double* sigmamap, int npix2, double*
arr<double> arrmap2(sigmamap, npix);
Healpix_Map <double> smap(arrmap2, RING);
gaussfit (hmap, smap, tabcat, ns, ncol, fwhm_arcmin, nfwhm , latcut, maxxi2 );
gaussfit (hmap, smap, tabcat, ns, ncol, fwhm_arcmin, nfwhm , latcut, maxxi2 , nofit);
arrmap = hmap.Map();
arrmap2.dealloc();
......
......@@ -38,5 +38,5 @@ void _make_grid_zone(double *map, int npix , int nside_in);
void _map2cov(double* map, int nin, int npix, double* mask, int npix2, double* cov, int ncross, int nz, double* nmode, int nz2, int reg=0);
void _combine_zone(double* map, int nin, int npix, double* mask, int npix2, double* cmap, int npix3, double* W, int nin2, int nz);
void _psdetection(double* map, int npix, double** cat, int* ncat, double fwhm_arcmin, double nsigma, int nside, double resolution, char* centers);
void _psgaussfit (double* cmap, int npix3, double* smap, int npix2, double* cat, int ns, int ncol,double fwhm_arcmin,double nfwhm ,double latcut,double maxxi2 ) ;
void _psgaussfit (double* cmap, int npix3, double* smap, int npix2, double* cat, int ns, int ncol,double fwhm_arcmin,double nfwhm ,double latcut,double maxxi2 , int nofit) ;
#endif /* !_SPHERELIB_MAP_H_ */
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment