Changeset 3788 in Sophya for trunk/Cosmo

Timestamp:

Jun 25, 2010, 7:35:17 PM (15 years ago)

Author:

ansari

Message:

Ajout classe de soustraction d'avant plans, Reza 25/06/2010

Location:

trunk/Cosmo/RadioBeam

Files:

• README (modified) (1 diff)
• applobe.cc (modified) (3 diffs)
• calcpk2.cc (modified) (8 diffs)
• cubedef.h (modified) (1 diff)
• fgndsub.cc (added)
• fgndsub.h (added)
• lobe.cc (modified) (2 diffs)
• lobe.h (modified) (2 diffs)
• makefile (modified) (2 diffs)
• mdish.cc (modified) (1 diff)
• mdish.h (modified) (2 diffs)

trunk/Cosmo/RadioBeam/README

@@ -25 +25 @@
 
 A.3/ List of other common files :
+  - fgndsub.h fgndsub.h : radio source, synchrotron subtraction class 
+  - lobe.h lobe.cc : 2D (k_alpha, k_delta) beam effect on 3D sky cubes
   - specpk.h specpk.cc : 3D power spectrum computation classes  
-  - lobe.h lobe.cc : 2D (k_alpha, k_delta) beam effect on 3D sky cubes
   - mdish.h mdish.cc : Classes for representing Multi-Dish interferometer configurations 
   - radutil.h radutil.cc : utilitaire de conversion a 21 cm

trunk/Cosmo/RadioBeam/applobe.cc

@@ -68 +68 @@
       pin >> incube;
     }
-    cout << incube;
+    incube.Show(); 
 
     double dxdeg = ThetaSizeDegre/(double)NTheta;
@@ -85 +85 @@
 
     cout << "applobe[2]: creating Four2DResponse and BeamEffect objects..." << endl;
-    // Dish de 64 m de diametre
-    double Diametre = 64.;
     H21Conversions conv;
     conv.setFrequency(Freq0MHz);
     double lambda = conv.getLambda();
-    Four2DResponse fresp(2, Diametre/lambda, Diametre/lambda);
+    Four2DResponse fresp(2, InterfArrayDiametre/lambda, InterfArrayDiametre/lambda);
     BeamEffect beam(fresp);
 
@@ -96 +94 @@
     beam.ApplyLobe3D(incube, dx, dy, Freq0MHz, dfreq);
 
-    cout << "applobe[4]: calling ReSample(incube, 0.25, 0.25, 1.) ... " << endl;
-    TArray< r_4 > outcube = beam.ReSample(incube, 0.25, 0.25, 1.);
+    cout << "applobe[4]: calling ReSample(incube, 0.5, 0.5, 1.) ... " << endl;
+    TArray< r_4 > outcube = beam.ReSample(incube, 0.5, 0.5, 1.);
     
-    cout << outcube;
+    outcube.Show(); 
     MeanSigma(outcube, mean, sigma);
     cout << " OutCube 3D- : Mean=" << mean << " Sigma=" << sigma <<  endl;

trunk/Cosmo/RadioBeam/calcpk2.cc

@@ -7 +7 @@
     R. Ansari , C. Magneville - Juin 2010 
 
-Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync InMapRadioSource convFacRsc OutPk [PixNoiseLevel]
+Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync InMapRadioSource convFacRsc OutPk 
+               [PixNoiseLevel] [TargetBeamArcmin] [NSigSrcThr]
 ---------------------------------------------------------------  */
 
@@ -20 +21 @@
 #include "specpk.h"
 #include "histats.h"
+#include "vector3d.h"
 
 #include "qhist.h" 
 #include "lobe.h" 
 #include "cubedef.h" 
+#include "fgndsub.h" 
+#include "radutil.h" 
 
 #include "histinit.h"
@@ -32 +36 @@
 
 typedef ThSDR48RandGen RandomGenerator ;
-//--- Declaration des fonctions 
-TArray<r_4> CleanForeground(TArray<r_4>& maplss, TArray<r_4>& mapsync, double freq0, double dfreq, 
-                            TArray<r_8>& synctemp, TArray<r_8>& specidx);
 
 //-------------------------------------------------------------------------
@@ -43 +44 @@
 {
   if (narg<6) {
-    cout << " Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync  OutPk [PixNoiseLevel] " << endl;
+    cout << " Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync  OutPk \n" 
+         << "        [PixNoiseLevel] [TargetBeamArcmin] [NSigSrcThr]" << endl;
     return 1;
   }
@@ -61 +63 @@
       fgaddnoise=true;
     }
+    
+    bool fgcorrbeam=true;
+    double tbeamarcmin=30.;
+    if (narg>7) {
+      tbeamarcmin=atof(arg[7]);
+      if (tbeamarcmin<1.e-6)  fgcorrbeam=false;
+    }
+    bool fgclnsrc=true;
+    double nsigsrc=5.;
+    if (narg>8) {
+      nsigsrc=atof(arg[8]);
+      if (nsigsrc<1.e-6)  fgclnsrc=false;
+    }
 
-
-    TArray<r_4> inmaplss, inmapsync;
-    const char * tits[2]={"LSS", "Sync"};
+    TArray<r_4> maplss, mapsync;
+    const char * tits[2]={"LSS", "Sync/RadioSrc"};
     for(int ks=0; ks<2; ks++)  {
       string& ppfname=inppflss;
       r_4 rfac=rfaclss;
-      TArray<r_4>* inmap=&inmaplss; 
-      if (ks==1) {  ppfname=inppfsync;  rfac=rfacsync;  inmap=&inmapsync; }
+      TArray<r_4>* inmap=&maplss; 
+      if (ks==1) {  ppfname=inppfsync;  rfac=rfacsync;  inmap=&mapsync; }
       cout << "calcpk2[" << ks+1 << "] : reading 3D map " << tits[ks] << " from file " << ppfname 
            << "  RenormFactor=" << rfac << endl;
@@ -82 +96 @@
     }
 
+    bool smo; 
+    if (!maplss.CompareSizes(mapsync,smo) ) {
+      cout << " calcpk2/ERROR sizes " << endl;
+      maplss.Show();  mapsync.Show(); 
+      return 99;
+    }
+    
+    TArray<r_4> skycube(mapsync);
+    skycube += maplss;
+
     if (fgaddnoise) { 
-      cout << " calcpk2: adding noise to LSS input cube ... " << endl;
-      BeamEffect::AddNoise(inmaplss, pixsignoise);
+      cout << " calcpk2: adding noise to skycube cube ... " << endl;
+      BeamEffect::AddNoise(skycube, pixsignoise);
     }
 
-    tm.Split(" After read ");
-    TArray<r_8> synctemp, specidx; 
-    cout << "calcpk2[3] : calling CleanForeground(...) " << endl;
-    double freq0 = Freq0MHz;
-    double dfreq = FreqSizeMHz/(double)NFreq;
-    TArray<r_4> inmap = CleanForeground(inmaplss, inmapsync, freq0, dfreq, synctemp, specidx);
     double mean, sigma;
-    MeanSigma(inmap, mean, sigma);
-    cout << " After cleaning: Mean=" << mean << " Sigma=" << sigma << endl;      
+    MeanSigma(skycube, mean, sigma);
+    cout << " input sky cube : Mean=" << mean << " Sigma=" << sigma << endl;      
+    tm.Split(" After input ");
+
+    H21Conversions conv;
+    conv.setFrequency(Freq0MHz);
+    double lambda = conv.getLambda();
+    Four2DResponse arep(2, InterfArrayDiametre/lambda, InterfArrayDiametre/lambda);
+    double DoL = 1.22/ArcminToRadian(tbeamarcmin);
+    Four2DResponse tbeam(2, DoL, DoL );
+    ForegroundCleaner  cleaner(arep, tbeam, skycube); 
+    if (fgcorrbeam) {
+      cout << "calcpk2[3.a] : calling cleaner.BeamCorrections() for target beam (" << tbeamarcmin << " arcmin)" 
+           << " Diam/Lambda=" << DoL << endl;
+      cleaner.BeamCorrections();
+    }
+    if (fgclnsrc) {
+      cout << "calcpk2[3.b] : calling cleaner.CleanPointSources() with threshold NSigma=" << nsigsrc  << endl;
+      cleaner.CleanPointSources(nsigsrc);
+    }
+
+    cout << "calcpk2[4] : calling cleaner.extractLSSCube(...) " << endl;
+    TArray<r_4> synctemp, specidx; 
+    TArray<r_4> exlss = cleaner.extractLSSCube(synctemp, specidx);
+
+    MeanSigma(exlss, mean, sigma);
+    cout << " After cleaning/extractLSS: Mean=" << mean << " Sigma=" << sigma << endl;      
     tm.Split(" After CleanForeground");
 
-    cout << "calcpk2[4] : computing 3D Fourier coefficients ... " << endl;
+    cout << "calcpk2[5] : computing 3D Fourier coefficients ... " << endl;
     FFTWServer ffts;
     TArray< complex<r_4> > four3d;
-    ffts.FFTForward(inmap, four3d);
+    ffts.FFTForward(exlss, four3d);
     tm.Split(" After FFTForward ");
 
-    cout << "calcpk2[5] : computing power spectrum ... " << endl;
+    cout << "calcpk2[6] : computing power spectrum ... " << endl;
     RandomGenerator rg;
     Four3DPk pkc(four3d, rg);
@@ -112 +155 @@
     tm.Split(" Done ComputePk ");    
     
-    cout << "calcpk2[4] : writing profile P(k) and foreground maps to  " << outname << endl;
+    cout << "calcpk2[7.a] : writing profile P(k)  to  " << outname << endl;
     POutPersist po(outname);
-    po << PPFNameTag("Pk") << hp;
-    po << PPFNameTag("Tsync") << synctemp;
-    po << PPFNameTag("async") << specidx;
+    po << hp;
+    outname = "fgm_" + outname;
+    cout << "calcpk2[7.b] : writing foreground maps to  " << outname << endl;
+    POutPersist pom(outname);
+    pom << PPFNameTag("Tsync") << synctemp;
+    pom << PPFNameTag("async") << specidx;
 
   }  // End of try bloc 
@@ -132 +178 @@
     rc = 97;
   }
-  cout << " ==== End of calcpk.cc program  Rc= " << rc << endl;
+  cout << " ==== End of calcpk2.cc program  Rc= " << rc << endl;
   return rc;    
 }
 
 
-/* --Fonction-- */
-TArray<r_4> CleanForeground(TArray<r_4>& maplss, TArray<r_4>& mapsync, double freq0, double dfreq, 
-                            TArray<r_8>& synctemp, TArray<r_8>& specidx)
-// Inputs : maplss, mapsyc, freq0, dfreq
-// Outputs : synctemp, specidx  (reconstructed foreground temperature and spectral index
-// Return_Array : foreground subtracted LSS signal 
-{
-  bool smo;
-  if (!maplss.CompareSizes(mapsync,smo) ) {
-    cout << " CleanForeground/ERROR sizes " << endl;
-    maplss.Show();  mapsync.Show(); 
-    throw ParmError("CleanForeground- maplss , mapsync not the same size !");
-  }
-  sa_size_t sz[5];   sz[0]=maplss.SizeX();  sz[1]=maplss.SizeY();
-  synctemp.SetSize(2, sz); 
-  specidx.SetSize(2, sz);
-  TArray<r_4>& omap=maplss;
-  Vector vlnf(maplss.SizeZ());
-  int nprt = 0;
-  // double freq0 : Frequence premier index en k (MHz)
-  // double dfreq :   // largeur en frequence de chaque plan (Mhz)  
-  for(sa_size_t i=0; i<maplss.SizeX(); i++) 
-    for(sa_size_t j=0; j<maplss.SizeY(); j++)  {
-      r_8 s1, sx, sx2, sy, sxy;
-      s1=sx=sx2=sy=sxy=0.;
-      for(sa_size_t k=0; k<maplss.SizeZ(); k++)  {
-        double lnf=log((double)k*dfreq+freq0);
-        vlnf(k)=lnf;
-        double ttot=(r_8)(mapsync(i,j,k)+maplss(i,j,k));
-        if (ttot < 1.e-5) continue;
-        double lntt=log(ttot);
-        s1+=1.;  sx+=lnf;  sx2+=(lnf*lnf);
-        sy+=lntt;    sxy+=(lnf*lntt);
-      }
-      double beta = (sx*sxy-sx2*sy)/(sx*sx-s1*sx2);
-      double alpha = (s1*sxy-sx*sy)/(s1*sx2-sx*sx);
-      double T0 = exp(beta+alpha*vlnf(0));
-      if ((i%6==0)&&(j%7==0)) 
-        cout << "CleanForeground[" << i << "," << j << "]: T0=" << T0 << " alpha=" << alpha 
-             << " (mapsync=" << mapsync(i,j,0) << " ... " << mapsync(i,j,125) << ")" << endl; 
-      synctemp(i,j) = T0;
-      specidx(i,j) = alpha;
-      for(sa_size_t k=0; k<maplss.SizeZ(); k++) {
-        r_4 fittedtemp = (r_4)(exp(beta+alpha*vlnf(k)));
-        omap(i,j,k) = mapsync(i,j,k)+maplss(i,j,k)-fittedtemp;
-      }
-    }
-  return omap;
-}

trunk/Cosmo/RadioBeam/cubedef.h

@@ -32 +32 @@
 static double sigPLidxSrc = 0.15;  // Sigma de la variation (gaussienne) de l'index  
 
+// Parametres pour la reponse en Fourier de l'instrument :
+static double InterfArrayDiametre = 64.;
+
 /*
 static sa_size_t NTheta=256;

trunk/Cosmo/RadioBeam/lobe.cc

@@ -17 +17 @@
 }
 
-/* --Methode-- */
-void BeamEffect::ApplyLobeK2D(TArray< complex<TF> >& fourAmp, double dkx, double dky, double lambda)
-//  dx, dy en radians, lambda en metres 
-{
-  fresp_.setLambda(lambda);
-  double kxx, kyy;
-  for(sa_size_t ky=0; ky<fourAmp.SizeY(); ky++) {
-    kyy =  (ky>fourAmp.SizeY()/2) ? -(double)(fourAmp.SizeY()-ky)*dky : (double)ky*dky; 
-    for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) {
-      kxx=(double)kx*dkx;
-      fourAmp(kx, ky) *= complex<TF>(fresp_(kxx, kyy), 0.);     
-    }
-  }
-  return;
-}
 
 
@@ -52 +37 @@
     ffts.FFTForward(slice, fourAmp);
     conv.setFrequency(f0+kz*df);
-    ApplyLobeK2D(fourAmp, dkx, dky, conv.getLambda());
+    fresp_.setLambda(conv.getLambda());
+    ApplyLobeK2D(fresp_, fourAmp, dkx, dky);
     ffts.FFTBackward(fourAmp, slice, true);
-    if (kz%10==0)  cout << "BeamEffect::ApplyLobe3D() done kz=" << kz << " / a.SizeZ()=" << a.SizeZ() << endl; 
+    if (kz%20==0)  cout << "BeamEffect::ApplyLobe3D() done kz=" << kz << " / a.SizeZ()=" << a.SizeZ() << endl; 
+  }
+  return;
+}
+
+/* --Methode-- */
+void BeamEffect::Correct2RefLobe(Four2DResponse& rep, TArray< TF >& a, double dx, double dy, double f0, double df)
+// dx, dy en radioans, f0, df en MHz
+{
+  Timer tm("BeamEffect::Correct2RefLobe");
+  FFTWServer ffts(true);                     
+  ffts.setNormalize(true); 
+  
+  H21Conversions conv;
+  
+  TArray< complex<TF> > fourAmp;
+  double dkx = DeuxPI/(double)a.SizeX()/dx;
+  double dky = DeuxPI/(double)a.SizeY()/dy;
+  
+  for(sa_size_t kz=0; kz<a.SizeZ(); kz++) {
+    TArray< TF > slice( a(Range::all(), Range::all(), kz) );
+    ffts.FFTForward(slice, fourAmp);
+    conv.setFrequency(f0+kz*df);
+    fresp_.setLambda(conv.getLambda());
+    Four2DRespRatio rratio(fresp_, rep);
+    ApplyLobeK2D(rratio, fourAmp, dkx, dky);
+    ffts.FFTBackward(fourAmp, slice, true);
+    if (kz%20==0)  cout << "BeamEffect::Correct2RefLobe() done kz=" << kz << " / a.SizeZ()=" << a.SizeZ() << endl; 
+  }
+  return;
+}
+
+/* --Methode-- */
+void BeamEffect::ApplyLobeK2D(Four2DResponse& rep, TArray< complex<TF> >& fourAmp, double dkx, double dky)
+//  dx, dy en radians, lambda en metres 
+{
+  double kxx, kyy;
+  for(sa_size_t ky=0; ky<fourAmp.SizeY(); ky++) {
+    kyy =  (ky>fourAmp.SizeY()/2) ? -(double)(fourAmp.SizeY()-ky)*dky : (double)ky*dky; 
+    for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) {
+      kxx=(double)kx*dkx;
+      fourAmp(kx, ky) *= complex<TF>(rep(kxx, kyy), 0.);        
+    }
   }
   return;

trunk/Cosmo/RadioBeam/lobe.h

@@ -1 +1 @@
 /*  ------------------------ Projet BAORadio -------------------- 
     Calcul de l'effet de lobe sur Carte2D (angleX,angleY) et 
-    cube 3D (angleX,angleY)
+    cube 3D (angleX,angleY, frquences)
     R. Ansari , C. Magneville - Juin 2010 
 ---------------------------------------------------------------  */
@@ -25 +25 @@
 class BeamEffect {
 public:
+  // Definition de l'objet avec la reponse en frequence de l'instrument 
   BeamEffect(Four2DResponse& resp);
-  void ApplyLobeK2D(TArray< complex<TF> >& f2d, double dkx, double dky, double lambda=1.);
+  // Applique l'effet de lobe au cube 3D (2 angles, frequence), pour chaque plan de frequence successivement  
   void ApplyLobe3D(TArray< TF >& a, double dx, double dy, double f0, double df);  
 
+  // Corrige de l'effet de l'effet de lobe, pour chaque plan de frequence, pour tout ramener au lobe defini par "rep" 
+  void Correct2RefLobe(Four2DResponse& rep, TArray< TF >& a, double dx, double dy, double f0, double df);
+
+  // Applique l'effet de lobe "rep" dans le plan de Fourier  pour une frequence (longueur d'onde) fixee
+  static void ApplyLobeK2D(Four2DResponse& rep, TArray< complex<TF> >& f2d, double dkx, double dky);
+
+  // Re-echntillonnage du cube 3D en appliquant les facteurs xfac,yfac,zfac selon chaque direction 
+  // fac = 2 ---> on double le nombre d'echantillon , fac=0.5 : un echantillon sur deux 
   static TArray< TF > ReSample(TArray< TF >& a, double xfac, double yfac, double zfac);
+
+  // On ajoute du bruit gaussien au cube 3D - espace des positions 
   static void AddNoise(TArray< TF >& a, double pixsignoise, bool fgcmsig=true);
 

trunk/Cosmo/RadioBeam/makefile

@@ -11 +11 @@
         rm Objs/* 
 
-PKGOLIST =  Objs/lobe.o Objs/specpk.o Objs/mdish.o Objs/qhist.o Objs/radutil.o 
-PKGHLIST =  lobe.h specpk.h mdish.h qhist.h radutil.h cubedef.h 
+PKGOLIST =  Objs/fgndsub.o Objs/lobe.o Objs/specpk.o Objs/mdish.o Objs/qhist.o Objs/radutil.o 
+PKGHLIST =  fgndsub.h lobe.h specpk.h mdish.h qhist.h radutil.h cubedef.h 
 
 ### les executables 
@@ -87 +87 @@
 
 ### les classes / fonctions 
-Objs/specpk.o : specpk.cc specpk.h  mdish.h qhist.h 
+Objs/fgndsub.o : fgndsub.cc $(PKGHLIST)
+        $(CXXCOMPILE) -o Objs/fgndsub.o fgndsub.cc
+
+Objs/specpk.o : specpk.cc $(PKGHLIST)
         $(CXXCOMPILE) -o Objs/specpk.o specpk.cc
 
-Objs/mdish.o : mdish.cc mdish.h qhist.h 
+Objs/mdish.o : mdish.cc $(PKGHLIST)
         $(CXXCOMPILE) -o Objs/mdish.o mdish.cc
 
-Objs/lobe.o : lobe.cc lobe.h mdish.h qhist.h 
+Objs/lobe.o : lobe.cc $(PKGHLIST)
         $(CXXCOMPILE) -o Objs/lobe.o lobe.cc
 

trunk/Cosmo/RadioBeam/mdish.cc

@@ -75 +75 @@
 }
 
+//---------------------------------------------------------------
+// -- Four2DRespRatio : rapport de la reponse entre deux objets Four2DResponse
+//---------------------------------------------------------------
+Four2DRespRatio::Four2DRespRatio(Four2DResponse& a, Four2DResponse& b)
+  : Four2DResponse(0, a.D(), a.D()), a_(a), b_(b) 
+{
+}
+
+double Four2DRespRatio::Value(double kx, double ky)
+{
+  double ra = a_.Value(kx,ky);
+  double rb = b_.Value(kx,ky);
+  if (rb<1.e-19) rb = 1.e-19;
+  return (ra/rb);
+}
+
+//---------------------------------------------------------------
 //--- Classe simple pour le calcul de rotation 
 class Rotation {

trunk/Cosmo/RadioBeam/mdish.h

@@ -55 +55 @@
 };
 
+
 // -- Four2DRespTable : Reponse tabulee instrumentale ds le plan k_x,k_y (angles theta,phi) 
 class Four2DRespTable : public  Four2DResponse {
@@ -63 +64 @@
   virtual double Value(double kx, double ky); 
   Histo2D hrep_;
+};
+
+
+// -- Four2DRespRatio: Retourne le rapport de la reponse entre deux objets Four2DResponse
+class Four2DRespRatio : public  Four2DResponse {
+public:
+  Four2DRespRatio(Four2DResponse& a, Four2DResponse& b);
+  // Return the ratio a.Value(kx,ky) / b.Value(kx, ky) - with protection against divide by zero 
+  virtual double Value(double kx, double ky); 
+  Four2DResponse& a_;
+  Four2DResponse& b_;
 };