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

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

File:

• trunk/Cosmo/RadioBeam/calcpk2.cc (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