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

Timestamp:

May 9, 2010, 12:47:58 AM (15 years ago)

Author:

cmv

Message:

speed-up FITS file read, cmv 09/05/2010

File:

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

trunk/Cosmo/SimLSS/cmvrvloscor.cc

@@ -11 +11 @@
 #include "timing.h"
 #include "dvlist.h"
-#include "ntuple.h"
+#include "histos.h"
 #include "fabtcolread.h"
 
@@ -27 +27 @@
 {
  int nthread = 1;
- SophyaInit();
- InitTim();
+ int_8 nhfill = 100000;
 
  if(narg<=2) {usage(); return -1;}
@@ -35 +34 @@
  try {
  //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
+
+ SophyaInit();
+ InitTim();
 
  cout<<"> read rho: "<<arg[1]<<endl;
@@ -52 +54 @@
  cout<<"Zref="<<Zref<<" Href="<<Href<<endl;
 
+ double dmin = min(Dx,min(Dy,Dz));
+ double nmax = max(Nx,max(Ny,Nz));
+ cout<<"dmin="<<dmin<<" nmax="<<nmax<<endl;
+ Histo hmpc(-dmin*nmax,dmin*nmax,4.*nmax);
+ nhfill = (int_8)Nx*(int_8)Ny*(int_8)Nz/nhfill; if(nhfill<=0) nhfill = 1;
+
  POutPersist pos("cmvrvloscor.ppf");
 
@@ -58 +66 @@
  TArray<GEN3D_TYPE>& rgen = fluct3d.GetRealArray();
  rgen = 0.;
- TVector<r_8> R(Nz), Rdis(Nz);
+ TVector<GEN3D_TYPE> R(Nz), V(Nz);
+ TVector<r_8> Rdis(Nz);
 
  cout<<"> filling redshift distorted cube"<<endl;
+ int_8 nread = 0;
  for(int i=0;i<Nx;i++) {
    if(i%(Nx/10)==0) cout<<"i="<<i<<endl;
  for(int j=0;j<Ny;j++) {
-   for(int l=0;l<Nz;l++) R(l) = f3dr.Read(l,j,i);
+   //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);
    Rdis = 0.;
    for(int l=0;l<Nz;l++) {
-     double d = (1.+Zref) / Href * f3dv.Read(l,j,i);
+     double d = (1.+Zref) / Href * V(l);
+     if(nread%nhfill==0) hmpc.Add(d);
      double lpd = (double)l + d/Dz; // valeur du deplacee
      // on repartit proportionnelement au recouvrement sur 2 pixels
@@ -80 +94 @@
  }
  PrtTim(">>>> End filling redshift distorted cube");
+ pos.PutObject(hmpc,"hmpc");
 
  fluct3d.ReComputeFourier();
@@ -131 +146 @@
 openppf cmvrvloscor.ppf
 
+disp hmpc
+
 n/plot hpkrec.val%x x>0 ! "nsta cpts logx"
 
@@ -144 +161 @@
 
 # proj selon kT (black), selon kZ (red)
-n/plot hpkrec2.val%sqrt(x*x+y*y) ! ! "nsta cpts"
+n/plot hpkrec2.val%sqrt(x*x+y*y) ! ! "nsta crossmarker3 logx"
  */