Changeset 3947 in Sophya for trunk/Cosmo/RadioBeam/pknoise.cc

Timestamp:

Feb 14, 2011, 12:58:29 AM (15 years ago)

Author:

ansari

Message:

Amelioration et modifs diverses lors de la redaction du papier, Reza 13/02/2011

File:

• trunk/Cosmo/RadioBeam/pknoise.cc (modified) (4 diffs)

trunk/Cosmo/RadioBeam/pknoise.cc

@@ -44 +44 @@
        << "    -renmax MaxValue (default : Do NOT renormalize 2D response value \n"   
        << "    -scut SCutValue (default=100.) \n"
-       << "    -ngen NGen (default=1) number of noise fourier amp generations \n"
+       << "    -ngen NGen (default=0) number of noise fourier amp generations \n"
+       << "       NGen=0 -> Call ComputeNoisePk(), else generate Fourier Amplitudes (random) \n"
        << "    -z redshift (default=0.7) \n"
        << "    -prt PrtLev,PrtModulo (default=0,10) " << endl;
@@ -74 +75 @@
   bool fgrenorm=false;
   double rmax=1.;
-  int NMAX = 1;
+  int NMAX = 0;
   double SCut=0.;
   double z_Redshift=0.7 ;  // 21 cm at z=0.7 -> 0.357 m  
@@ -141 +142 @@
       cout << "pknoise[1]: initializing Four2DRespTable from file" << resptblname << endl;
       resptbl.readFromPPF(resptblname);
+      cout << "pknoise[1.b] Four2DRespTable.LambdaRef=" << resptbl.getLambdaRef() << " setLambda(" 
+           << lambda << ")" << endl;
+      resptbl.setLambda(lambda);
       arep_p=&resptbl;
       if (fgrenorm) {
-        cout << " pknoise[1.b] call to resptbl.renormalize(" << rmax << ")"; 
+        cout << "pknoise[1.c] call to resptbl.renormalize(" << rmax << ")"; 
         double omax=resptbl.renormalize(rmax);
         cout << " ... Old Max=" << omax << endl;
@@ -159 +163 @@
 
     cout << " pknoise[3]: Computing Noise P(k) using PkNoiseCalculator ..." << endl;
-    PkNoiseCalculator pkn(m3d, *(arep_p), SCut, NMAX, tits.c_str());
-    pkn.SetPrtLevel(prtlev,prtmod);
-    HProf hpn = pkn.Compute();
-    po << hpn;
- 
+    if (NMAX>0) {
+      PkNoiseCalculator pkn(m3d, *(arep_p), SCut, NMAX, tits.c_str());
+      pkn.SetPrtLevel(prtlev,prtmod);
+      HProf hpn = pkn.Compute();
+      po << hpn;
+    }
+    else {
+      Histo fracmodok;
+      DataTable dtnoise;
+      HProf hpn = m3d.ComputeNoisePk(*(arep_p),fracmodok,dtnoise,SCut);
+      po << hpn;
+      string outfile2 = "x"+outfile;
+      POutPersist po2(outfile2);
+      HProf h1dnoise=arep_p->GetProjNoiseLevel();
+      HProf h1drep=arep_p->GetProjResponse();
+      po2 << PPFNameTag("dtnoise") << dtnoise;
+      po2 << PPFNameTag("hpnoise") << hpn;
+      po2 << PPFNameTag("fracmodok") << fracmodok;
+      po2 << PPFNameTag("h1dnoise") << h1dnoise;
+      po2 << PPFNameTag("h1drep") << h1drep;
+    }
     rc = 0;
   }  // End of try bloc