source: Sophya/trunk/Cosmo/SimLSS/cmvrvloscor.cc@ 3802

Last change on this file since 3802 was 3802, checked in by cmv, 15 years ago

minor cosmetics, cmv 22/07/2010
File size: 11.8 KB
Line 
1#include "sopnamsp.h"
2#include "machdefs.h"
3#include <iostream>
4#include <stdlib.h>
5#include <stdio.h>
6#include <string.h>
7#include <math.h>
8#include <unistd.h>
9
10#include "sophyainit.h"
11#include "timing.h"
12#include "dvlist.h"
13#include "histos.h"
14#include "fabtcolread.h"
15#include "fftwserver.h"
16
17#include "constcosmo.h"
18#include "geneutils.h"
19#include "genefluct3d.h"
20// set simul = 6_0p0_780
21// cmvrvloscorf -o rvlos_${simul}.ppf -K 75 -S ginit3d_${simul}_r.fits ginit3d_${simul}_rv.fits
22// cmvrvloscorf -o rvlos_${simul}.ppf -n 1,30 -K 75 -S -2 ginit3d_${simul}_r.fits ginit3d_${simul}_rv.fits
23// cmvrvloscorf -o rvlos_${simul}.ppf -n 1,30 -2 ginit3d_${simul}_r.fits ginit3d_${simul}_rv.fits
24
25void usage(void);
26void usage(void)
27{
28cout
29<<"cmvrvloscor [options] rho.fits vlos.fits"<<endl
30<<"-n nplany,nhfill: process one Y plane every \"nplany\" (def:1(all))"<<endl
31<<" fill histos with \"nhfill\" los (def:25)"<<endl
32<<"-K npix: compute correlation R*V at +/- npix pixels (def: no)"<<endl
33<<" (very time consuming!!!)"<<endl
34<<"-S: compute cross-power spectrum of V*conj(R) (def: no)"<<endl
35<<"-N: do not create 3D cube and recompute 1D and 2D spectra (def: no do-it !)"<<endl
36<<"-p: do not distribute dRho/Rho in pixel, just take the nearest (def: no, distribute rho)"<<endl
37<<"-2: compute 2D projection fpr dRho and dRho(corrected) (def=no)"<<endl
38<<"-o: output ppf file name (def=\"cmvrvloscor.ppf\")"<<endl
39<<endl;
40}
41
42int main(int narg,char *arg[])
43{
44 int nthread = 1, nplany=1, nhfilllos = 25, npixcor = 0;
45 bool docube=true, dopk = false, do2d = false, dodistrib=true;
46 string fnppf = "cmvrvloscor.ppf";
47
48 // --- Decodage des arguments
49 char c;
50 while((c = getopt(narg,arg,"hn:K:SNp2o:")) != -1) {
51 switch (c) {
52 case 'n' :
53 sscanf(optarg,"%d,%d",&nplany,&nhfilllos);
54 if(nplany<=0) nplany = 1;
55 if(nhfilllos<=0) nhfilllos = 0;
56 break;
57 case 'K' :
58 npixcor = atoi(optarg);
59 break;
60 case 'S' :
61 dopk = true;
62 break;
63 case 'N' :
64 docube = false;
65 break;
66 case '2' :
67 do2d = true;
68 break;
69 case 'p' :
70 dodistrib = false;
71 break;
72 case 'o' :
73 fnppf = optarg;
74 break;
75 case 'h' :
76 default :
77 usage(); return -1;
78 }
79 }
80 if(optind>=narg-1) {usage(); return -1;}
81
82 //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
83 try {
84 //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
85
86 SophyaInit();
87 InitTim();
88
89 // --- open FITS files (dRho/Rho and Vlos)
90 cout<<"> read rho: "<<arg[optind]<<endl;
91 FitsImg3DRead f3dr(arg[optind],0,5);
92 cout<<"> read vlos: "<<arg[optind+1]<<endl;
93 FitsImg3DRead f3dv(arg[optind+1],0,5);
94 long Nx = f3dr.ReadKeyL("Nx");
95 long Ny = f3dr.ReadKeyL("Ny");;
96 long Nz = f3dr.ReadKeyL("Nz");;
97 cout<<"N: x="<<Nx<<" y="<<Ny<<" z="<<Nz<<endl;
98 double Dx = f3dr.ReadKey("Dx");
99 double Dy = f3dr.ReadKey("Dy");
100 double Dz = f3dr.ReadKey("Dz");
101 cout<<"D: x="<<Dx<<" y="<<Dy<<" z="<<Dz<<endl;
102 double Zref = f3dr.ReadKey("ZREF");
103 double Href = f3dr.ReadKey("HREF");
104 cout<<"Zref="<<Zref<<" Href="<<Href<<endl;
105
106 double dmin = min(Dx,min(Dy,Dz));
107 long nmax = max(Nx,max(Ny,Nz));
108 cout<<"dmin="<<dmin<<" nmax="<<nmax<<endl;
109 Histo hmpc(-dmin*nmax/10.,dmin*nmax/10.,nmax);
110
111 POutPersist pos(fnppf.c_str());
112 DVList dvlcor;
113
114 // --- Create a Cube for analysis
115 GeneFluct3D *fluct3d = NULL;
116 TArray<GEN3D_TYPE>* rgen = NULL;
117 if(docube) {
118 cout<<"...Create and fill 3D cube"<<endl;
119 fluct3d = new GeneFluct3D(Nx,Ny,Nz,Dx,Dy,Dz,nthread,2);
120 fluct3d->Print();
121 rgen = &(fluct3d->GetRealArray());
122 *rgen = 0.;
123 cout<<"rgen: size [1]="<<rgen->SizeX()
124 <<" [2]="<<rgen->SizeY()
125 <<" [3]="<<rgen->SizeZ()<<endl;
126 cout<<"pkgen: size [1]="<<fluct3d->GetComplexArray().SizeX()
127 <<" [2]="<<fluct3d->GetComplexArray().SizeY()
128 <<" [3]="<<fluct3d->GetComplexArray().SizeZ()<<endl;
129 }
130
131 // --- Vector for real-space correlation computation
132 int imil = Nz-1;
133 dvlcor("imil") = (int_4)imil;
134 TVector<int_4> nKsi;
135 TVector<r_8> Ksirv, Ksirvc, Ksirr, Ksirrc;
136 if(npixcor>0) {
137 Ksirv.ReSize(2*Nz-1); Ksirv = 0.;
138 Ksirvc.ReSize(2*Nz-1); Ksirvc = 0.;
139 Ksirr.ReSize(2*Nz-1); Ksirr = 0.;
140 Ksirrc.ReSize(2*Nz-1); Ksirrc = 0.;
141 nKsi.ReSize(2*Nz-1); nKsi = 0;
142 cout<<"...Compute R*V correlation on +/-"<<npixcor<<" px"<<endl;
143 }
144
145 // --- Vector for PK cross-correlation computation
146 int npk = 0;
147 TVector< complex<r_4> > FR, FV;
148 TVector< complex<r_8> > pkvr, FRdis;
149 TVector<r_8> pkr, pkrc;
150 FFTWServer fftserv;
151 if(dopk) cout<<"...compute V*conj(R) cross-power spectrum"<<endl;
152
153 // --- Read and process data
154 TVector<r_4> R(Nz), V(Nz);
155 TVector<r_8> Rdis(Nz);
156 if(nplany>Ny) nplany = Ny;
157 cout<<"...Will read one Y plane every "<<nplany<<endl;
158 if(nhfilllos) {
159 cout<<"...Fill Mpc displacement histo with "<<nhfilllos<<" los"<<endl;
160 nhfilllos = int((double)Nx*Ny/nplany/nhfilllos + 0.5);
161 if(nhfilllos<=0) nhfilllos = 1;
162 cout<<" -> fill one los every "<<nhfilllos<<endl;
163 }
164
165 cout<<">>> filling redshift distorted cube"<<endl;
166 int nlosread = 0;
167 for(int i=0;i<Nx;i++) {
168 if(i%(Nx/10)==0) cout<<" i="<<i<<endl;
169 TMatrix<r_4> M2d, M2dv, M2dc;
170 if(do2d && (i==0 || i==Nx/2 || i==Nx-1)) {
171 M2d.ReSize(Ny,Nz); M2d = 0.;
172 M2dv.ReSize(Ny,Nz); M2dv = 0.;
173 M2dc.ReSize(Ny,Nz); M2dc = 0.;
174 }
175 for(int j=0;j<Ny;j+=nplany) {
176 bool fhis = false;
177 if(nhfilllos) if(nlosread%nhfilllos==0) fhis = true;
178 //for(int l=0;l<Nz;l++) R(l) = f3dr.Read(l,j,i);
179 //for(int l=0;l<Nz;l++) V(l) = f3dv.Read(l,j,i);
180 f3dr.Read(j,i,R);
181 f3dv.Read(j,i,V);
182 Rdis = 0.;
183 // Calcul du champ R redshift distordu
184 for(int l=0;l<Nz;l++) {
185 double d = V(l) / Href;
186 if(fhis) hmpc.Add(d);
187 double lpd = (double)l + d/Dz; // valeur du deplacee
188 if(dodistrib) {
189 // on repartit proportionellement au recouvrement sur 2 pixels
190 long l1 = long(lpd); // pixel de gauche
191 long l2 = l1 + 1; // pixel de droite
192 lpd -= (double)l1; // recouvrement du pixel du dessus
193 if(l1>=0 && l1<Nz) Rdis(l1) += R(l) * (1.-lpd);
194 if(l2>=0 && l2<Nz) Rdis(l2) += R(l) * lpd;
195 } else { // take nearest pixel
196 long l1 = long(lpd + 0.5);
197 if(l1>=0 && l1<Nz) Rdis(l1) += R(l);
198 }
199 }
200 // On remplit eventuellement les matrices 2D
201 if(do2d && M2d.Size()>0)
202 for(int l=0;l<Nz;l++) {M2d(j,l) = R(l); M2dv(j,l) = V(l); M2dc(j,l) = Rdis(l);}
203 // On remplit le cube avec le champ R redshift distordu
204 if(fluct3d) for(int l=0;l<Nz;l++) (*rgen)(l,j,i) += Rdis(l);
205 // Calcul eventuel de la fonction de correlation R*V
206 if(npixcor>0) {
207 for(long l1=0;l1<Nz;l1++) {
208 for(long l2=max(0L,l1-npixcor);l2<min(Nz,l1+npixcor);l2++) {
209 int lc = imil+(l2-l1);
210 Ksirr(lc) += R(l1)*R(l2);
211 Ksirrc(lc) += Rdis(l1)*R(l2);
212 Ksirv(lc) += R(l1)*V(l2);
213 Ksirvc(lc) += Rdis(l1)*V(l2);
214 nKsi(lc)++;
215 }
216 }
217 }
218 // Cross-power spectrum computation
219 if(dopk) {
220 fftserv.FFTForward(V,FV);
221 int nf = FV.Size();
222 if(pkvr.Size()<=0) {
223 cout<<"...Create vector for cross-power spectrum computation"<<endl;
224 pkvr.ReSize(nf); pkvr = complex<r_8>(0.);
225 pkr.ReSize(nf); pkr = 0.;
226 pkrc.ReSize(nf); pkrc = 0.;
227 }
228 fftserv.FFTForward(R,FR);
229 for(int l=0;l<nf;l++) {
230 pkvr(l) += FV(l)*conj(FR(l));
231 pkr(l) += norm(FR(l));
232 }
233 fftserv.FFTForward(Rdis,FRdis);
234 for(int l=0;l<nf;l++) pkrc(l) += norm(FRdis(l));
235 npk++;
236 }
237 nlosread++;
238 }
239 if(do2d && M2d.Size()>0) {
240 char str[64];
241 sprintf(str,"mx_%d",i);
242 pos.PutObject(M2d,str);
243 sprintf(str,"mxv_%d",i);
244 pos.PutObject(M2dv,str);
245 sprintf(str,"mxc_%d",i);
246 pos.PutObject(M2dc,str);
247 }
248 }
249
250 cout<<"Number of processed los: "<<nlosread<<" / "<<Nx*Ny<<endl;
251 dvlcor("nlosread") = (int_4)nlosread;
252 if(hmpc.NEntries()>0) {
253 hmpc.Show();
254 pos.PutObject(hmpc,"hmpc");
255 }
256 if(Ksirr.Size()>0) {
257 for(int l=0;l<Ksirr.Size();l++) if(nKsi(l)>0) {
258 Ksirr(l) /= nKsi(l);
259 Ksirrc(l) /= nKsi(l);
260 Ksirv(l) /= nKsi(l);
261 Ksirvc(l) /= nKsi(l);
262 }
263 pos.PutObject(Ksirr,"ksirr");
264 pos.PutObject(Ksirrc,"ksirrc");
265 pos.PutObject(Ksirv,"ksirv");
266 pos.PutObject(Ksirvc,"ksirvc");
267 pos.PutObject(nKsi,"nksi");
268 }
269 if(npk>0) {
270 pkvr /= (double)npk;
271 pkr /= (double)npk;
272 pkrc /= (double)npk;
273 pos.PutObject(pkvr,"pkvr");
274 pos.PutObject(pkr,"pkr");
275 pos.PutObject(pkrc,"pkrc");
276 }
277 PrtTim(">>>> End filling redshift distorted cube");
278
279 // --- Fourier transform 3D cube and compute 1D and 2D power spectra
280 if(fluct3d) {
281 cout<<">>> Fourier transform 3D cube and compute 1D and 2D power spectra"<<endl;
282 // do the FFT for spectrum analysis
283 fluct3d->ReComputeFourier();
284 PrtTim(">>>> End ReComputing spectrum");
285
286 // Compute 1D spectrum
287 cout<<endl<<"\n--- Computing final 1D spectrum"<<endl;
288 double dkmin = fluct3d->GetKincMin(), knyqmax = fluct3d->GetKmax();
289 long nherr = long(knyqmax/dkmin+0.5);
290 cout<<"\nFor HistoErr: d="<<dkmin<<" max="<<knyqmax<<" n="<<nherr<<endl;
291 HistoErr hpkrec(0.,knyqmax,nherr); hpkrec.Zero();
292 hpkrec.ReCenterBin(); hpkrec.Show();
293 fluct3d->ComputeSpectrum(hpkrec);
294 pos.PutObject(hpkrec,"hpkrec");
295 PrtTim(">>>> End Computing final spectrum");
296
297 // Compute 2D spectrum
298 cout<<"\n--- Computing final 2D spectrum"<<endl;
299 double dktmin = fluct3d->GetKTincMin(), ktnyqmax = fluct3d->GetKTmax();
300 double dkzmin = fluct3d->GetKinc()[2], kznyqmax = fluct3d->GetKnyq()[2];
301 long nherrt = long(ktnyqmax/dktmin+0.5), nherrz = long(kznyqmax/dkzmin+0.5);
302 cout<<"For Histo2DErr: d="<<dktmin<<","<<dkzmin
303 <<" max="<<ktnyqmax<<","<<kznyqmax<<" n="<<nherrt<<","<<nherrz<<endl;
304 Histo2DErr hpkrec2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
305 hpkrec2.ReCenterBin(); hpkrec2.Zero(); hpkrec2.Show();
306 fluct3d->ComputeSpectrum2D(hpkrec2);
307 pos.PutObject(hpkrec2,"hpkrec2");
308 PrtTim(">>>> End Computing final 2D spectrum");
309 }
310
311 // --- end of job, write objects in ppf
312 pos.PutObject(dvlcor,"dvlcor");
313 if(fluct3d) delete fluct3d;
314
315 //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
316 } catch (PException& exc) {
317 cerr<<"cmvrvloscor.cc catched PException"<<exc.Msg()<<endl;
318 return 77;
319 } catch (std::exception& sex) {
320 cerr << "cmvrvloscor.cc std::exception :"
321 << (string)typeid(sex).name() << "\n msg= "
322 << sex.what() << endl;
323 return 78;
324 } catch (...) {
325 cerr << "cmvrvloscor.cc catched unknown (...) exception " << endl;
326 return 79;
327 }
328 //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
329
330 return 0;
331}
332
333/*
334openppf cmvrvloscor.ppf
335
336disp hmpc
337
338# cross-correlation
339disp nksi
340set imil ${dvlcor.imil}
341n/plot ksirv.val%n-${imil} ! ! "nsta cpts"
342n/plot ksirvc.val%n-${imil} ! ! "nsta cpts same red"
343
344n/plot ksirr.val%n-${imil} ! ! "nsta cpts"
345n/plot ksirrc.val%n-${imil} ! ! "nsta cpts same red"
346
347# cross-power spectrum
348n/plot pkr.val%n ! ! "nsta cpts logx"
349n/plot pkrc.val%n ! ! "nsta cpts logx same red"
350
351n/plot pkvr.val%n ! ! "nsta cpts logx"
352
353# reconstructed 1D power spectrum
354n/plot hpkrec.val%x x>0 ! "nsta cpts logx"
355
356# reconstructed 2D power spectrum
357imag hpkrec2
358addoval 0 0 0.05 0.05 "green" false
359addoval 0 0 0.1 0.1 "green" false
360addoval 0 0 0.25 0.25 "green" false
361addoval 0 0 0.5 0.5 "green" false
362x = ${hpkrec2.xmax} / 2.
363addoval 0 0 $x $x "green" false
364x = ${hpkrec2.ymax} / 2.
365addoval 0 0 $x $x "green" false
366
367# proj selon kT (black), selon kZ (red)
368n/plot hpkrec2.val%sqrt(x*x+y*y) ! ! "nsta crossmarker3 logx"
369
370# les matrices 2D
371set n 0
372disp mx_$n
373newwin
374disp mxc_$n
375
376defscript tom
377del m2
378c++exec TMatrix<r_8> m2(hpkrec2.NBinX(),hpkrec2.NBinY()); \
379for(int i=0;i<hpkrec2.NBinX();i++) \
380for(int j=0;j<hpkrec2.NBinY();j++) m2(i,j) = hpkrec2(i,j); \
381KeepObj(m2); cout<<"OK"<<endl;
382endscript
383
384###### Comparaison avec le spectre initial
385cd /home
386mkdir ok
387mkdir gene
388cd /home
389
390cd /ok
391openppf cmvrvloscor.ppf
392tom
393cd /home
394
395cd /gene
396openppf ginit3d_?_?p?_???.ppf
397tom
398cd /home
399
400h/oper / /ok/m2 /gene/m2 mdiff
401disp mdiff
402 */
Note: See TracBrowser for help on using the repository browser.