Changeset 3821 in Sophya for trunk/Cosmo

Timestamp:

Jul 30, 2010, 12:35:57 PM (15 years ago)

Author:

cmv

Message:

type d'interpolation pour les spectres tabules, amelioration de cmvrvloscor.cc, cmv 30/07/2010

Location:

trunk/Cosmo/SimLSS

Files:

• cmvginit3d.cc (modified) (1 diff)
• cmvrvloscor.cc (modified) (10 diffs)
• cmvtstpk.cc (modified) (1 diff)

trunk/Cosmo/SimLSS/cmvginit3d.cc

@@ -250 +250 @@
    pkzt->ReadCAMB(fn,h100tab,ztab);
    pkzt->SetGrowthFactor(&growth);
+   pkzt->SetInterpTyp(2);
    pkzt->SetZ(zref);
    pkz = pkzt;

trunk/Cosmo/SimLSS/cmvrvloscor.cc

@@ -19 +19 @@
 #include "genefluct3d.h"
 // set simul = 6_0p0_780
-// 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
+// cmvrvloscorf -o rvlos_${simul}.ppf -n 1,1 -z -1,0 -H 50 -2 -S 10000 -s 1. ginit3d_${simul}_r.fits ginit3d_${simul}_rv.fits
 
 void usage(void);
@@ -28 +26 @@
 cout
 <<"cmvrvloscor [options] rho.fits vlos.fits"<<endl
-<<"-n nplany,nhfill: process one Y plane every \"nplany\" (def:1(all))"<<endl
-<<"                  fill histos with \"nhfill\" los (def:25)"<<endl
-<<"-K npix: compute correlation R*V at +/- npix pixels (def: no)"<<endl
-<<"         (very time consuming!!!)"<<endl
-<<"-S: compute cross-power spectrum of V*conj(R) (def: no)"<<endl
+<<"-n nplany,nplanx: process one Y (X) plane every \"nplany\" (\"nplanx\") (def:1,1(all))"<<endl
+<<"-z nplanz,iplanz0: process only \"nplanz\" Z plane beginnig at plane \"iplanz0\" (def:all Z planes)"<<endl
+<<"-H nhfill: fill displacement histo with \"nhfill\" los (def:25)"<<endl
+<<"-s scaldis: scale the displacement by factor \"scaldis\" (def: 1.)"<<endl
+<<"-p: do not distribute dRho/Rho in pixel, just take the nearest (def: no, distribute rho)"<<endl
+<<"-S npkmax: compute cross-power spectrum of V*conj(R) with \"npkmax\" los (def: no)"<<endl
+<<"-2: compute 2D projection for dRho and dRho(corrected) for 3 Y-Z plans (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
+<<"-o ppfout.ppf: output ppf file name (def=\"cmvrvloscor.ppf\")"<<endl
 <<endl;
 }
@@ -42 +40 @@
 int main(int narg,char *arg[])
 {
- int nthread = 1, nplany=1, nhfilllos = 25, npixcor = 0;
- bool docube=true, dopk = false, do2d = false, dodistrib=true;
+ int nthread = 1, nplany=1, nplanx=1, nplanz=-1, iplanz0=0, nhfilllos = 25, npkmax = 0;
+ double scaldis = 1.;
+ bool docube=true, do2d = false, dodistrib=true;
  string fnppf = "cmvrvloscor.ppf";
  
  // --- Decodage des arguments
  char c;
- while((c = getopt(narg,arg,"hn:K:SNp2o:")) != -1) {
+ while((c = getopt(narg,arg,"hn:H:S:Np2o:s:z:")) != -1) {
   switch (c) {
   case 'n' :
-    sscanf(optarg,"%d,%d",&nplany,&nhfilllos);
+    sscanf(optarg,"%d,%d",&nplany,&nplanx);
     if(nplany<=0) nplany = 1;
+    if(nplanx<=0) nplanx = 1;
+    break;
+  case 'z' :
+    sscanf(optarg,"%d,%d",&nplanz,&iplanz0);
+    break;
+  case 'H' :
+    nhfilllos = atoi(optarg);
     if(nhfilllos<=0) nhfilllos = 0;
     break;
-  case 'K' :
-    npixcor = atoi(optarg);
-    break;
   case 'S' :
-    dopk = true;
+    npkmax = atoi(optarg);
     break;
   case 'N' :
     docube = false;
     break;
+  case 'p' :
+    dodistrib = false;
+    break;
   case '2' :
     do2d = true;
     break;
-  case 'p' :
-    dodistrib = false;
-    break;
   case 'o' :
     fnppf = optarg;
+    break;
+  case 's' :
+    scaldis = atof(optarg);
     break;
   case 'h' :
@@ -104 +110 @@
  cout<<"Zref="<<Zref<<" Href="<<Href<<endl;
 
- double dmin = min(Dx,min(Dy,Dz));
- long nmax = max(Nx,max(Ny,Nz));
- cout<<"dmin="<<dmin<<" nmax="<<nmax<<endl;
- Histo hmpc(-dmin*nmax/10.,dmin*nmax/10.,nmax);
-
  POutPersist pos(fnppf.c_str());
- DVList dvlcor;
+ DVList dvllos;
+
+ // --- Read and process data
+ if(nplany>Ny || nplany<=0) nplany = Ny;
+ if(nplanx>Nx || nplanx<=0) nplanx = Nx;
+ cout<<"...Will read one Y plane every "<<nplany<<endl;
+ cout<<"             one X plane every "<<nplanx<<endl;
+ int nlospred = Nx*Ny/nplany/nplanx;
+ cout<<"   predicted number of los to be read "<<nlospred<<endl;
+ if(nplanz>Nz || nplanz<=0) nplanz = Nz;
+ if(iplanz0<0)iplanz0 = 0;
+ if(iplanz0>=Nz) {cout<<"Error: iplanz0="<<iplanz0<<" out of range"<<endl; return -2;}
+ if(iplanz0+nplanz>Nz) nplanz = Nz - iplanz0;
+ cout<<"   process "<<nplanz<<" Z planes in ["<<iplanz0<<","<<iplanz0+nplanz<<"["<<endl;
+ TVector<r_4> V4read(Nz);
+ TVector<r_4> R(nplanz), V(nplanz);
+ TVector<r_8> Rdis(nplanz);
+
+ // --- Histo des deplacements
+ Histo* hmpc = NULL;
+ if(nhfilllos>0) {
+   cout<<"...Fill Mpc displacement histo with "<<nhfilllos<<" los"<<endl;
+   double dmin = min(Dx,min(Dy,Dz)); if(fabs(scaldis)>0.) dmin *= fabs(scaldis);
+   long nmax = max(Nx,max(Ny,Nz));
+   hmpc = new Histo(-dmin*nmax/10.,dmin*nmax/10.,nmax);
+   cout<<"   dmin="<<dmin<<" nmax="<<nmax<<endl;
+   nhfilllos = int((double)nlospred/nhfilllos + 0.5);
+   if(nhfilllos<=0) nhfilllos = 1;
+   cout<<"   -> fill one los every "<<nhfilllos<<endl;
+ }
+
+ // --- Vector for PK cross-correlation computation
+ if(npkmax<0) npkmax = Nx*Ny;
+ TVector< complex<r_4> > FR, FV;
+ TVector< complex<r_8> > pkvr, FRdis;
+ TVector<r_8> pkr, pkrc;
+ FFTWServer fftserv;
+ fftserv.setNormalize(false); // pour accord avec GeneFluct3D
+ if(npkmax>0) {
+   cout<<"...compute R*conj(R) V*conj(V) V*conj(R) cross-power spectrum with "<<npkmax<<" los"<<endl;
+   npkmax = int((double)nlospred/npkmax + 0.5);
+   if(npkmax<=0) npkmax = 1;
+   cout<<"   -> fill one los every "<<npkmax<<endl;
+ }
+ int npkav = 0;
+ double dkzpk = 2.*M_PI/(nplanz*Dz);  // Mpc^-1
 
  // --- Create a Cube for analysis
@@ -117 +163 @@
  if(docube) {
    cout<<"...Create and fill 3D cube"<<endl;
-   fluct3d = new GeneFluct3D(Nx,Ny,Nz,Dx,Dy,Dz,nthread,2);
+   fluct3d = new GeneFluct3D(Nx,Ny,nplanz,Dx,Dy,Dz,nthread,2);
    fluct3d->Print();
    rgen = &(fluct3d->GetRealArray());
@@ -129 +175 @@
  }
 
- // --- Vector for real-space correlation computation
- int imil = Nz-1;
- dvlcor("imil") = (int_4)imil;
- TVector<int_4> nKsi;
- TVector<r_8> Ksirv, Ksirvc, Ksirr, Ksirrc;
- if(npixcor>0) {
-   Ksirv.ReSize(2*Nz-1);  Ksirv  = 0.;
-   Ksirvc.ReSize(2*Nz-1); Ksirvc = 0.;
-   Ksirr.ReSize(2*Nz-1);  Ksirr  = 0.;
-   Ksirrc.ReSize(2*Nz-1); Ksirrc = 0.;
-   nKsi.ReSize(2*Nz-1);   nKsi   = 0;
-   cout<<"...Compute R*V correlation on +/-"<<npixcor<<" px"<<endl;
- }
-
- // --- Vector for PK cross-correlation computation
- int npk = 0;
- TVector< complex<r_4> > FR, FV;
- TVector< complex<r_8> > pkvr, FRdis;
- TVector<r_8> pkr, pkrc;
- FFTWServer fftserv;
- if(dopk) cout<<"...compute V*conj(R) cross-power spectrum"<<endl;
-
- // --- Read and process data
- TVector<r_4> R(Nz), V(Nz);
- TVector<r_8> Rdis(Nz);
- if(nplany>Ny) nplany = Ny;
- cout<<"...Will read one Y plane every "<<nplany<<endl;
- if(nhfilllos) {
-   cout<<"...Fill Mpc displacement histo with "<<nhfilllos<<" los"<<endl;
-   nhfilllos = int((double)Nx*Ny/nplany/nhfilllos + 0.5);
-   if(nhfilllos<=0) nhfilllos = 1;
-   cout<<"   -> fill one los every "<<nhfilllos<<endl;
- }
-
- cout<<">>> filling redshift distorted cube"<<endl;
+ // --- Do the jobs !!!
+ cout<<">>> filling redshift distorted cube, scale displacement by "<<scaldis<<endl;
  int nlosread = 0;
- for(int i=0;i<Nx;i++) {
-   if(i%(Nx/10)==0) cout<<"   i="<<i<<endl;
+ int lpri = nlospred/25; if(lpri==0) lpri = 1;
+ for(int i=0;i<Nx;i+=nplanx) {
    TMatrix<r_4> M2d, M2dv, M2dc;
    if(do2d && (i==0 || i==Nx/2 || i==Nx-1)) {
-     M2d.ReSize(Ny,Nz); M2d = 0.;
-     M2dv.ReSize(Ny,Nz); M2dv = 0.;
-     M2dc.ReSize(Ny,Nz); M2dc = 0.;
+     M2d.ReSize(Ny,nplanz); M2d = 0.;
+     M2dv.ReSize(Ny,nplanz); M2dv = 0.;
+     M2dc.ReSize(Ny,nplanz); M2dc = 0.;
    }
  for(int j=0;j<Ny;j+=nplany) {
-   bool fhis = false;
-   if(nhfilllos) if(nlosread%nhfilllos==0) fhis = true;
-   //for(int l=0;l<Nz;l++) R(l) = f3dr.Read(l,j,i);
-   //for(int l=0;l<Nz;l++) V(l) = f3dv.Read(l,j,i);
-   f3dr.Read(j,i,R);
-   f3dv.Read(j,i,V);
+   if(nlosread%lpri==0) cout<<"   nlosread="<<nlosread<<" (i="<<i<<",j="<<j<<")"<<endl;
+   bool fhis = (hmpc != NULL && nlosread%nhfilllos==0) ? true: false;
+   //for(int l=0;l<nplanz;l++) R(l) = f3dr.Read(iplanz0+l,j,i);
+   //for(int l=0;l<nplanz;l++) V(l) = f3dv.Read(iplanz0+l,j,i);
+   f3dr.Read(j,i,V4read);
+     for(int l=0;l<nplanz;l++) R(l) = V4read(iplanz0+l);
+   f3dv.Read(j,i,V4read);
+     for(int l=0;l<nplanz;l++) V(l) = V4read(iplanz0+l);
    Rdis = 0.;
    // Calcul du champ R redshift distordu
-   for(int l=0;l<Nz;l++) {
-     double d = V(l) / Href;
-     if(fhis) hmpc.Add(d);
+   for(int l=0;l<nplanz;l++) {
+     double d = scaldis * V(l) / Href;
+     if(fhis) hmpc->Add(d);
      double lpd = (double)l + d/Dz; // valeur du deplacee
      if(dodistrib) { 
@@ -191 +206 @@
        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(l1>=0 && l1<nplanz) Rdis(l1) += R(l) * (1.-lpd);
+       if(l2>=0 && l2<nplanz) Rdis(l2) += R(l) * lpd;
      } else { // take nearest pixel
        long l1 = long(lpd + 0.5);
-       if(l1>=0 && l1<Nz) Rdis(l1) += R(l);
+       if(l1>=0 && l1<nplanz) Rdis(l1) += R(l);
      }
    }
    // On remplit eventuellement les matrices 2D
    if(do2d && M2d.Size()>0)
-     for(int l=0;l<Nz;l++) {M2d(j,l) = R(l); M2dv(j,l) = V(l); M2dc(j,l) = Rdis(l);}
+     for(int l=0;l<nplanz;l++) {M2d(j,l) = R(l); M2dv(j,l) = V(l); M2dc(j,l) = Rdis(l);}
    // On remplit le cube avec le champ R redshift distordu
-   if(fluct3d) for(int l=0;l<Nz;l++) (*rgen)(l,j,i) += Rdis(l);
-   // Calcul eventuel de la fonction de correlation R*V
-   if(npixcor>0) {
-     for(long l1=0;l1<Nz;l1++) {
-       for(long l2=max(0L,l1-npixcor);l2<min(Nz,l1+npixcor);l2++) {
-         int lc = imil+(l2-l1);
-         Ksirr(lc)  += R(l1)*R(l2);
-         Ksirrc(lc) += Rdis(l1)*R(l2);
-         Ksirv(lc)  += R(l1)*V(l2);
-         Ksirvc(lc) += Rdis(l1)*V(l2);
-         nKsi(lc)++;
-       }
-     }
-   }
+   if(fluct3d) for(int l=0;l<nplanz;l++) (*rgen)(l,j,i) += Rdis(l);
    // Cross-power spectrum computation
-   if(dopk) {
+   if(npkmax>0 && nlosread%npkmax==0) {
      fftserv.FFTForward(V,FV);
      int nf = FV.Size();
      if(pkvr.Size()<=0) {
-       cout<<"...Create vector for cross-power spectrum computation"<<endl;
-       pkvr.ReSize(nf); pkvr = complex<r_8>(0.);
+       cout<<"...Create vector for cross-power spectrum computation sz="<<nf<<endl;
        pkr.ReSize(nf);  pkr = 0.;
        pkrc.ReSize(nf); pkrc = 0.;
+       pkvr.ReSize(nf); pkvr = complex<r_8>(0.);
      }
      fftserv.FFTForward(R,FR);
+     fftserv.FFTForward(Rdis,FRdis);
      for(int l=0;l<nf;l++) {
+       pkr(l)  += norm(FR(l));
+       pkrc(l) += norm(FRdis(l));
        pkvr(l) += FV(l)*conj(FR(l));
-       pkr(l)  += norm(FR(l));
      }
-     fftserv.FFTForward(Rdis,FRdis);
-     for(int l=0;l<nf;l++) pkrc(l) += norm(FRdis(l));
-     npk++;
+     npkav++;
    }
    nlosread++;
- }
+ }   // fin boucle sur ny
    if(do2d && M2d.Size()>0) {
      char str[64];
@@ -246 +248 @@
      pos.PutObject(M2dc,str);
    }
- }
+ }   // fin boucle sur nx
 
  cout<<"Number of processed los: "<<nlosread<<" / "<<Nx*Ny<<endl;
- dvlcor("nlosread") = (int_4)nlosread;
- if(hmpc.NEntries()>0) {
-   hmpc.Show();
-   pos.PutObject(hmpc,"hmpc");
- }
- if(Ksirr.Size()>0) {
-   for(int l=0;l<Ksirr.Size();l++) if(nKsi(l)>0) {
-       Ksirr(l)  /= nKsi(l);
-       Ksirrc(l) /= nKsi(l);
-       Ksirv(l)  /= nKsi(l);
-       Ksirvc(l) /= nKsi(l);
-     }
-   pos.PutObject(Ksirr,"ksirr");
-   pos.PutObject(Ksirrc,"ksirrc");
-   pos.PutObject(Ksirv,"ksirv");
-   pos.PutObject(Ksirvc,"ksirvc");
-   pos.PutObject(nKsi,"nksi");
- }
- if(npk>0) {
-   pkvr /= (double)npk;
-   pkr  /= (double)npk;
-   pkrc /= (double)npk;
-   pos.PutObject(pkvr,"pkvr");
+ if(hmpc != NULL) {
+   hmpc->Show();
+   if(hmpc->NEntries()>0) pos.PutObject(*hmpc,"hmpc");
+ }
+ if(npkav>0) {
+   double nn = Dx*Dy*Dz/double(nplanz)/double(npkav); // moyenne + normalisation / hpkrec
+   cout<<"Number of averaged spectra: npkav="<<npkav<<" norm="<<nn<<endl;
+   pkr  *= nn;
+   pkrc *= nn;
+   pkvr *= nn;
    pos.PutObject(pkr,"pkr");
    pos.PutObject(pkrc,"pkrc");
+   pos.PutObject(pkvr,"pkvr");
  }
  PrtTim(">>>> End filling redshift distorted cube");
@@ -310 +300 @@
 
  // --- end of job, write objects in ppf
- pos.PutObject(dvlcor,"dvlcor");
- if(fluct3d) delete fluct3d;
+ dvllos("nlosread") = (int_4)nlosread;
+ dvllos("nlospred") = (int_4)nlospred;
+ dvllos("Nx") = (int_4)Nx;
+ dvllos("Ny") = (int_4)Ny;
+ dvllos("Nz") = (int_4)Nz;
+ dvllos("Dx") = (r_8)Dx;
+ dvllos("Dy") = (r_8)Dy;
+ dvllos("Dz") = (r_8)Dz;
+ dvllos("nplanz") = (int_4)nplanz;
+ dvllos("iplanz0") = (int_4)iplanz0;
+ dvllos("nplany") = (int_4)nplany;
+ dvllos("nplanx") = (int_4)nplanx;
+ dvllos("npkav") = (int_4)npkav;
+ dvllos("dkzpk") = (r_8)dkzpk;
+ pos.PutObject(dvllos,"dvllos");
+ if(fluct3d != NULL) delete fluct3d;
+ if(hmpc != NULL) delete hmpc;
 
  //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
@@ -334 +339 @@
 openppf cmvrvloscor.ppf
 
+c++exec dvllos.Print();
+set dkz ${dvllos.dkzpk}
+
 disp hmpc
 
-# cross-correlation
-disp nksi
-set imil ${dvlcor.imil}
-n/plot ksirv.val%n-${imil} ! ! "nsta cpts"
-n/plot ksirvc.val%n-${imil} ! ! "nsta cpts same red"
-
-n/plot ksirr.val%n-${imil} ! ! "nsta cpts"
-n/plot ksirrc.val%n-${imil} ! ! "nsta cpts same red"
-
 # cross-power spectrum
-n/plot pkr.val%n ! ! "nsta cpts logx"
-n/plot pkrc.val%n ! ! "nsta cpts logx same red"
-
-n/plot pkvr.val%n ! ! "nsta cpts logx"
+n/plot pkr.val%n*${dkz} n>0&&val>0 ! "nsta cpts logx logy"
+n/plot pkrc.val%n*${dkz} n>0&&val>0 ! "nsta cpts logx logy same red"
+
+n/plot pkvr.val%n*${dkz} n>0&&val>0 ! "nsta cpts logx logy"
 
 # reconstructed 1D power spectrum
-n/plot hpkrec.val%x x>0 ! "nsta cpts logx"
+n/plot hpkrec.val%x x>0&&val>0 ! "nsta cpts logx logy"
 
 # reconstructed 2D power spectrum

trunk/Cosmo/SimLSS/cmvtstpk.cc

@@ -114 +114 @@
    pkcamb = new PkTabulate;
    pkcamb->ReadCAMB(fcmbfile,h100,0.);
-   pkcamb->SetInterpTyp(1);
+   pkcamb->SetInterpTyp(2);
    if(pkcamb->NPoints()==0) {delete pkcamb; pkcamb = NULL;}
  }