Changeset 3788 in Sophya for trunk/Cosmo/RadioBeam/calcpk2.cc

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

File:

• trunk/Cosmo/RadioBeam/calcpk2.cc (modified) (8 diffs)

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;
-}