Changeset 3830 in Sophya

Timestamp:

Aug 4, 2010, 1:59:10 PM (15 years ago)

Author:

ansari

Message:

Amelioration soustraction des avant-plans par l'introduction d'un fit polynome deg 2 sur ln(Temp)=f(ln(freq)), Reza 04/08/2010

Location:

trunk/Cosmo/RadioBeam

Files:

• calcpk2.cc (modified) (5 diffs)
• fgndsub.cc (modified) (6 diffs)
• fgndsub.h (modified) (1 diff)
• subtractradsrc.cmd (modified) (5 diffs)

trunk/Cosmo/RadioBeam/calcpk2.cc

@@ -43 +43 @@
 int main(int narg, const char* arg[])
 {
-  if (narg<6) {
-    cout << " Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync  OutPkFile \n" 
-         << "        [PixNoiseLevel] [Diameter/Four2DRespTableFile] [TargetBeamArcmin] [NSigSrcThr]" << endl;
-    return 1;
+  if ( (narg<6)||((narg>1)&&(strcmp(arg[1],"-h")==0)) ) {
+    cout << " Usage: calcpk2 InMapLSS convFacLSS InMapFgnd convFacFgnd OutPkFile \n" 
+         << "        [PixNoiseLevel] [D_Dish/Four2DRespTableFile CorBeamDoL] \n"
+         << "        [NSigSrcThr] [P2/P1] [RecMapFile] " << endl;
+    if ((narg>1)&&(strcmp(arg[1],"-h")==0)) {
+      cout << "- InMapLSS: Input 3D LSS cube (PPF file name) \n "
+           << "- convFacLSS: LSS cube conversion factor to mK (milliKelvin) \n"
+           << "- InMapFgnd: Input 3D foreground cube (PPF file name) \n"
+           << "- convFacFgnd: Foreground cube conversion factor to mK (milliKelvin) \n"
+           << "- PixNoiseLevel: White noise level per pixel (mK) (default=0.) \n" 
+           << "- D_Dish/Four2DRespTableFile: Dish diameter or 2D (u,v) plane response (PPF file name) \n"
+           << "- CorBeamDoL: Beam correction target Diameter/Lambda \n" 
+           << "    These two parameters are used to correct for beam effect for a \n"
+           << "    target beam (independent of frequency) defined by D/Lambda \n"
+           << "    DoL = 100 --> beam ~ 35 arcmin (D=30m @ z~0.5) \n" 
+           << "    default : no beam correction applied \n " 
+           << " - NSigSrcThr: Point source cleaning, Nb_Sigmas on stacked 2D temperature \n" 
+           << "    default : no point source cleaning, use NSigSrcThr ~ 3..5 \n"
+           << " - P2/P1: 2nd/first degree polynomial fit on ln(Temp) = f(ln(freq)) \n "
+           << "    foreground subtraction. default is P2 \n" 
+           << "- RecMapFile: output PPF file for reconstructed foreground template \n" 
+           << "    (Temperature,SpectralIndex) and extracted LSS cube \n"
+         << endl;
+      return 1;
+    }
+    else cout << "   calcpk2 -h for detailed usage " << endl;
+    return 2;
   }
   Timer tm("calcpk2");
@@ -79 +102 @@
       }
     }
-    double tbeamDoL=135;
+    double tbeamDoL=0.;
     if (narg>8) {
       tbeamDoL=atof(arg[8]);
@@ -90 +113 @@
       if (nsigsrc<1.e-6)  fgclnsrc=false;
     }
-
+    bool fgpoly2=true;  // true -> soustraction polynome degre 2
+    if ((narg>0)&&(strcmp(arg[10],"P1")==0)) fgpoly2=false;      
+    bool fgsavemaps=false;
+    string outmap_ppfname="extlss.ppf";
+    if (narg>11) {
+      outmap_ppfname=arg[11];
+      fgsavemaps=true;
+    }
+    
     TArray<r_4> maplss, mapsync;
     const char * tits[2]={"LSS", "Sync/RadioSrc"};
@@ -165 +196 @@
     cout << "calcpk2[4] : calling cleaner.extractLSSCube(...) " << endl;
     TArray<r_4> synctemp, specidx; 
-    TArray<r_4> exlss = cleaner.extractLSSCube(synctemp, specidx);
+    TArray<r_4> exlss; 
+    if (fgpoly2) exlss = cleaner.extractLSSCubeP2(synctemp, specidx);
+    else  exlss = cleaner.extractLSSCubeP1(synctemp, specidx);
 
     MeanSigma(exlss, mean, sigma);
@@ -188 +221 @@
 
     tm.Split(" Done ComputePk ");    
-    
+    { 
     cout << "calcpk2[7.a] : writing profile P(k)  to  " << outname << endl;
     POutPersist po(outname);
     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;
+    }
+    if (fgsavemaps) {
+      cout << "calcpk2[7.b] : writing foreground maps and extracted LSS to  " << outmap_ppfname << endl;
+      POutPersist pom(outmap_ppfname);
+      pom << PPFNameTag("Tsync") << synctemp;
+      pom << PPFNameTag("async") << specidx;
+      pom << PPFNameTag("extlss") << exlss;
+    }
 
   }  // End of try bloc 

trunk/Cosmo/RadioBeam/fgndsub.cc

@@ -9 +9 @@
 #include "cubedef.h"
 #include "matharr.h"
+#include "poly.h"
 
 #include "ctimer.h"
@@ -90 +91 @@
 
 /* --Methode-- */
-TArray< TF >  ForegroundCleaner::extractLSSCube(TArray< TF >& synctemp, TArray< TF >& specidx)
-{
-  Timer tm("ForegroundCleaner::extractLSSCube");
+TArray< TF >  ForegroundCleaner::extractLSSCubeP1(TArray< TF >& synctemp, TArray< TF >& specidx)
+{
+  Timer tm("ForegroundCleaner::extractLSSCubeP1");
 // Inputs : maplss, mapsyc, freq0, dfreq
 // Outputs : synctemp, specidx  (reconstructed foreground temperature and spectral index
@@ -130 +131 @@
       bool fgnan = false;
       if (!isfinite(alpha)||(!isfinite(beta))) {
-        cout << "extractLSSCube[" << i << "," << j << "]/ Not finite alpha, beta - (mapsync="
+        cout << "extractLSSCubeP1[" << i << "," << j << "]/ Not finite alpha, beta - (mapsync="
              << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl; 
         alpha=beta=-999.;
@@ -140 +141 @@
         
         if ((axp1<-70.)||(axp1>70.)||(axp2<-70.)||(axp2>70.)) {
-        cout << "extractLSSCube[" << i << "," << j << "] BAD alpha=" << alpha 
+        cout << "extractLSSCubeP1[" << i << "," << j << "] BAD alpha=" << alpha 
              << " beta=" << beta << " T0=" << T0 << " - (mapsync="
              << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl; 
@@ -147 +148 @@
       }
       if ((i%imodprt==0)&&(j%jmodprt==0)) 
-        cout << "extractLSSCube[" << i << "," << j << "]: T0=" << T0 << " alpha=" << alpha 
+        cout << "extractLSSCubeP1[" << i << "," << j << "]: T0=" << T0 << " alpha=" << alpha 
              << " (mapsync=" << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl; 
       synctemp(i,j) = T0;
@@ -161 +162 @@
       }
     }
-  cout << " ForegroundCleaner::extractLSSCube() - NbNan alpha/beta=" << nbinfini << " NbBAD =" << nbbad << endl;
+  cout << " ForegroundCleaner::extractLSSCubeP1() - NbNan alpha/beta=" << nbinfini << " NbBAD =" << nbbad << endl;
   return omap;
 }
 
+/*
+static inline val_polyn2(double alpha, double beta, double gamma, double x)
+{
+  return (beta+alpha*x+gamma*x*x);
+}
+*/
+
+/* --Methode-- */
+TArray< TF >  ForegroundCleaner::extractLSSCubeP2(TArray< TF >& synctemp, TArray< TF >& specidx)
+{
+  Timer tm("ForegroundCleaner::extractLSSCubeP2");
+// Inputs : maplss, mapsyc, freq0, dfreq
+// Outputs : synctemp, specidx  (reconstructed foreground temperature and spectral index
+// Return_Array : foreground subtracted LSS signal 
+  sa_size_t sz[5];   sz[0]=skycube_.SizeX();  sz[1]=skycube_.SizeY();
+  synctemp.SetSize(2, sz); 
+  specidx.SetSize(2, sz);
+  TArray<r_4> omap; 
+  omap.SetSize(skycube_, true);
+  Vector vlnf(skycube_.SizeZ());
+  Vector vlnT(skycube_.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 k=0; k<skycube_.SizeZ(); k++)  
+    vlnf(k)=log((double)k*dfreq_+freq0_);
+  vlnf -= vlnf(0);
+  //  cout << " DBG*extractLSSCubeP2 vlnf(0)=" << vlnf(0) << " vlnf(1)=" << vlnf(1) 
+  //       << "vlnf(last)=" << vlnf(vlnf.Size()-1) << endl;
+  sa_size_t nbinfini=0;
+  sa_size_t nbbad=0;
+  sa_size_t imodprt=skycube_.SizeX()/6;
+  sa_size_t jmodprt=skycube_.SizeY()/6;
+  for(sa_size_t i=0; i<skycube_.SizeX(); i++) 
+    for(sa_size_t j=0; j<skycube_.SizeY(); j++)  {
+      vlnT = -12.;
+      Poly polyn;
+      for(sa_size_t k=0; k<skycube_.SizeZ(); k++)  {
+        double lnf = vlnf(k);
+        double ttot=(r_8)(skycube_(i,j,k));
+        if (ttot < 1.e-5) continue;
+        vlnT(k)=log(ttot);
+      }
+      polyn.Fit(vlnf,vlnT,2);
+      double beta = polyn[0];
+      double alpha = polyn[1];
+      double gamma = polyn[2];
+      double T0 = exp(polyn(vlnf(0))); // exp( val_polyn2(alpha, beta, gamma, vlnf(0)) ); 
+
+      bool fgnan = false;
+      if (!isfinite(alpha)||(!isfinite(beta))||(!isfinite(gamma))) {
+        cout << "extractLSSCubeP2[" << i << "," << j << "]/ Not finite alpha, beta - (mapsync="
+             << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl; 
+        alpha=beta=gamma=-999.;
+        fgnan = true;  nbinfini++;
+      }
+      else {
+        double axp1 = polyn(vlnf(0));
+        double axp2 = polyn(vlnf(vlnf.Size()-1));
+        
+        if ((axp1<-70.)||(axp1>70.)||(axp2<-70.)||(axp2>70.)) {
+          cout << "extractLSSCubeP2[" << i << "," << j << "] BAD alpha=" << alpha 
+            //         << " beta=" << beta << " gamma=" << gamma << " T0=" << T0 << " axp1=" << axp1 << " axp2=" << axp2 
+               << " beta=" << beta << " gamma=" << gamma << " T0=" << T0 << " - (mapsync=" << skycube_(i,j,0) 
+               << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl;
+          fgnan = true;  nbbad++;
+        }
+      }
+      if ((i%imodprt==0)&&(j%jmodprt==0)) 
+        cout << "extractLSSCubeP2[" << i << "," << j << "]: T0=" << T0 << " alpha=" << alpha << " gamma=" << gamma  
+             << " (mapsync=" << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl; 
+      synctemp(i,j) = T0;
+      specidx(i,j) = alpha;
+      if (fgnan) {
+        for(sa_size_t k=0; k<skycube_.SizeZ(); k++)   omap(i,j,k) = 0.;
+      }
+      else {
+        for(sa_size_t k=0; k<skycube_.SizeZ(); k++) {
+          r_4 fittedtemp = (r_4)( exp(polyn(vlnf(k))) );
+          omap(i,j,k) = skycube_(i,j,k)-fittedtemp;
+        }
+      }
+    }
+  cout << " ForegroundCleaner::extractLSSCubeP2() - NbNan alpha/beta=" << nbinfini << " NbBAD =" << nbbad << endl;
+  return omap;
+}
+

trunk/Cosmo/RadioBeam/fgndsub.h

@@ -29 +29 @@
   int CleanNegatives(TF seuil=1.e-6);
   int CleanPointSources(double nsigmas=5.);
-  TArray< TF > extractLSSCube(TArray< TF >& synctemp, TArray< TF >& specidx);
+  // Ajustement d'une droite (a x + b) sur ln(Temp) = f(ln(Freq)) 
+  TArray< TF > extractLSSCubeP1(TArray< TF >& synctemp, TArray< TF >& specidx);
+  // Ajustement d'une parabole (a x^2 + b x + c) sur ln(Temp) = f(ln(Freq)) 
+  TArray< TF > extractLSSCubeP2(TArray< TF >& synctemp, TArray< TF >& specidx);
 
   Four2DResponse& arrep_;   // Array/Instrument beam response

trunk/Cosmo/RadioBeam/subtractradsrc.cmd

@@ -84 +84 @@
 
 # 4.c/ Extract LSS P(k) from Foreground+LSS+noise , after cleaning/subtraction without beam 
-csh> ./Objs/calcpk2 lsscube.ppf 0.2 fgndcube.ppf 1000 subpk.ppf 0.35 50. 0. 0.
+csh> ./Objs/calcpk2 lsscube.ppf 0.2 fgndcube.ppf 1000 subpk.ppf 0.35 50. 0. 0. P2
 # 4.d / Extract LSS P(k) from Foreground+LSS+noise and beam effect, without beam correction  
-csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.2 fgndcube_lobe.ppf 1000 subpklobe.ppf 0.35 50. 0. 3. 
+csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.2 fgndcube_lobe.ppf 1000 subpklobe.ppf 0.35 50. 0. 3. P2 
 # 4.e / Extract LSS P(k) from Foreground+LSS+noise and beam effect - correcting to a beam of Diam/Lambda = 130 
-csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.2 fgndcube_lobe.ppf 1000 subpkcorlobe.ppf 0.35 50. 130. 3.
+csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.2 fgndcube_lobe.ppf 1000 subpkcorlobe.ppf 0.35 50. 130. 3. P2 
+#  Or using a linear fit for foreground subtraction (old version)
+csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.2 fgndcube_lobe.ppf 1000 subpkcorlobep1.ppf 0.35 50. 130. 3. P1
 # 4.f / Estimate residual noise from Foreground removal : 
-csh> ./Objs/calcpk2 zerolss.ppf 0. fgndcube_lobe.ppf 1000 subpknolss.ppf 0.35 50. 130. 3.
+csh> ./Objs/calcpk2 zerolss.ppf 0. fgndcube_lobe.ppf 1000 subpknolss.ppf 0.35 50. 130. 3. P2
 csh> ./Objs/calcpk2 zerolss.ppf 0. fgndcube_lobe.ppf 1000 subpknolssnocor.ppf 0.35 50. 0. 3.
 
@@ -103 +105 @@
 openppf lsspklobewn.ppf 
 openppf subpklobe.ppf 
+
 openppf subpkcorlobe.ppf 
+openppf subpknolss.ppf
 
-openppf subpknolss.ppf
-openppf subpknolssnocor.ppf
+# openppf subpknolssnocor.ppf
 
 disp lsspk 'logx logy nsta xylimits=0.005,2.,4e-11,8e-6 gold'
@@ -113 +116 @@
 settitle ' Pk[LSS] - without normalisation' ' ' 'font=helvetica,bold,16 black' 
 
-#  Calcul du volume total en Mpc^3
-set VOL 3*3*3*360*360*256
-plot2d lsspklobe x val*$VOL 1 'same nsta orange'
-plot2d lsspklobe x val*$VOL 1 'same nsta orange'
 
 disp fgndpk 'logx logy nsta xylimits=0.005,2.,1e-10,1. navyblue'
@@ -123 +122 @@
 disp lsspk 'same nsta gold'
 disp lsspklobewn 'same nsta siennared' 
-settitle 'Pk[LSS] , Pk[Foreground] and lobe effect (Dish D=50 m)' ' ' 'font=helvetica,bold,18'
+# settitle 'Pk[LSS] , Pk[Foreground] and lobe effect (Dish D=50 m)' ' ' 'font=helvetica,bold,18'
+settitle 'Pk[LSS] , Pk[Foreground=GSM] and lobe effect (Dish D=50 m)' ' ' 'font=helvetica,bold,18'
+
 set lines ( 'Pk[Foreground]'  'Pk[fgnd]*Lobe' 'Pk[fgnd]*Lobe/Corrected' 'Pk[LSS]' 'Pk[LSS]*Lobe+Noise' )
 set cols ( navyblue blue skyblue gold siennared ) 
 textdrawer lines cols 'font=helvetica,bold,16 frame'
+
 
 disp lsspk 'logx logy nsta xylimits=0.005,2.,4e-9,4e-5 gold'
@@ -141 +143 @@
 plot2d subpknolss x val*$VOL 1 'same nsta cpts marker=box,5 green' 
 
+# settitle 'Recovered Pk[LSS] In=LSS+(GSM) (D=50 m)' ' ' 'font=helvetica,bold,18'
 settitle 'Recovered Pk[LSS] In=LSS+(Haslam+North20cm) (D=50 m)' ' ' 'font=helvetica,bold,18'
 setaxelabels 'k (Mpc^-1)   h=0.7' 'P(k)    (mK^2 Mpc^3)' 'font=helvetica,bolditalic,16'