Changeset 3799 in Sophya for trunk/Cosmo/SimLSS/cmvrvloscor.cc

Timestamp:

Jun 30, 2010, 6:23:54 PM (15 years ago)

Author:

cmv

Message:

Vlos * by (1+z), Dlos(com) do not * by (1+z), cmv 30/06/2010

File:

• trunk/Cosmo/SimLSS/cmvrvloscor.cc (modified) (7 diffs)

trunk/Cosmo/SimLSS/cmvrvloscor.cc

@@ -19 +19 @@
 #include "genefluct3d.h"
 // set simul = 6_0p0_780
-// cmvrvloscorf -K 75 -S ginit3d_${simul}_r.fits ginit3d_${simul}_rv.fits
-// cmvrvloscorf -n 1,30 -K 75 -S -2 ginit3d_${simul}_r.fits ginit3d_${simul}_rv.fits
-// cmvrvloscorf -n 1,30 -2 ginit3d_${simul}_r.fits ginit3d_${simul}_rv.fits
+// cmvrvloscorf -o rvlos_${simul}.ppf -K 75 -S ginit3d_${simul}_r.fits ginit3d_${simul}_rv.fits
+// cmvrvloscorf -o rvlos_${simul}.ppf -n 1,30 -K 75 -S -2 ginit3d_${simul}_r.fits ginit3d_${simul}_rv.fits
+// cmvrvloscorf -o rvlos_${simul}.ppf -n 1,30 -2 ginit3d_${simul}_r.fits ginit3d_${simul}_rv.fits
 
 void usage(void);
@@ -34 +34 @@
 <<"-S: compute cross-power spectrum of V*conj(R) (def: no)"<<endl
 <<"-N: do not create 3D cube and recompute 1D and 2D spectra (def: no do-it !)"<<endl
+<<"-p: do not distribute dRho/Rho in pixel, just take the nearest (def: no, distribute rho)"<<endl
 <<"-2: compute 2D projection fpr dRho and dRho(corrected) (def=no)"<<endl
+<<"-o: output ppf file name (def=\"cmvrvloscor.ppf\")"<<endl
 <<endl;
 }
@@ -41 +43 @@
 {
  int nthread = 1, nplany=1, nhfilllos = 25, npixcor = 0;
- bool docube=true, dopk = false, do2d = false;
+ bool docube=true, dopk = false, do2d = false, dodistrib=true;
+ string fnppf = "cmvrvloscor.ppf";
  
  // --- Decodage des arguments
  char c;
- while((c = getopt(narg,arg,"hn:K:SN2")) != -1) {
+ while((c = getopt(narg,arg,"hn:K:SNp2o:")) != -1) {
   switch (c) {
   case 'n' :
@@ -63 +66 @@
   case '2' :
     do2d = true;
+    break;
+  case 'p' :
+    dodistrib = false;
+    break;
+  case 'o' :
+    fnppf = optarg;
     break;
   case 'h' :
@@ -100 +109 @@
  Histo hmpc(-dmin*nmax,dmin*nmax,4*nmax);
 
- POutPersist pos("cmvrvloscor.ppf");
+ POutPersist pos(fnppf.c_str());
  DVList dvlcor;
 
@@ -173 +182 @@
    // Calcul du champ R redshift distordu
    for(int l=0;l<Nz;l++) {
-     double d = (1.+Zref) / Href * V(l);
+     double d = V(l) / Href;
      if(fhis) hmpc.Add(d);
      double lpd = (double)l + d/Dz; // valeur du deplacee
-     // on repartit proportionellement au recouvrement sur 2 pixels
-     long l1 = long(lpd); // pixel de gauche
-     long l2 = l1 + 1;  // pixel de droite
-     lpd -= (double)l1;  // recouvrement du pixel du dessus
-     if(l1>=0 && l1<Nz) Rdis(l1) += R(l) * (1.-lpd);
-     if(l2>=0 && l2<Nz) Rdis(l2) += R(l) * lpd;
+     if(dodistrib) { 
+       // on repartit proportionellement au recouvrement sur 2 pixels
+       long l1 = long(lpd); // pixel de gauche
+       long l2 = l1 + 1;  // pixel de droite
+       lpd -= (double)l1;  // recouvrement du pixel du dessus
+       if(l1>=0 && l1<Nz) Rdis(l1) += R(l) * (1.-lpd);
+       if(l2>=0 && l2<Nz) Rdis(l2) += R(l) * lpd;
+     } else { // take nearest pixel
+       long l1 = long(lpd + 0.5);
+       if(l1>=0 && l1<Nz) Rdis(l1) += R(l);
+     }
    }
    // On remplit eventuellement les matrices 2D
@@ -351 +365 @@
 n/plot hpkrec2.val%sqrt(x*x+y*y) ! ! "nsta crossmarker3 logx"
 
+defscript tom
+del m2
+c++exec TMatrix<r_8> m2(hpkrec2.NBinX(),hpkrec2.NBinY()); \
+for(int i=0;i<hpkrec2.NBinX();i++) \
+for(int j=0;j<hpkrec2.NBinY();j++) m2(i,j) = hpkrec2(i,j); \
+KeepObj(m2); cout<<"OK"<<endl;
+endscript
+
 # les matrices 2D
 set n 0