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

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

Vlos * by (1+z), Dlos(com) do not * by (1+z), cmv 30/06/2010
File size: 11.5 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,dmin*nmax,4*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, M2dc;
170 if(do2d && (i==0 || i==Nx/2 || i==Nx-1)) {
171 M2d.ReSize(Ny,Nz); M2d = 0.;
172 M2dc.ReSize(Ny,Nz); M2dc = 0.;
173 }
174 for(int j=0;j<Ny;j+=nplany) {
175 bool fhis = false;
176 if(nhfilllos) if(nlosread%nhfilllos==0) fhis = true;
177 //for(int l=0;l<Nz;l++) R(l) = f3dr.Read(l,j,i);
178 //for(int l=0;l<Nz;l++) V(l) = f3dv.Read(l,j,i);
179 f3dr.Read(j,i,R);
180 f3dv.Read(j,i,V);
181 Rdis = 0.;
182 // Calcul du champ R redshift distordu
183 for(int l=0;l<Nz;l++) {
184 double d = V(l) / Href;
185 if(fhis) hmpc.Add(d);
186 double lpd = (double)l + d/Dz; // valeur du deplacee
187 if(dodistrib) {
188 // on repartit proportionellement au recouvrement sur 2 pixels
189 long l1 = long(lpd); // pixel de gauche
190 long l2 = l1 + 1; // pixel de droite
191 lpd -= (double)l1; // recouvrement du pixel du dessus
192 if(l1>=0 && l1<Nz) Rdis(l1) += R(l) * (1.-lpd);
193 if(l2>=0 && l2<Nz) Rdis(l2) += R(l) * lpd;
194 } else { // take nearest pixel
195 long l1 = long(lpd + 0.5);
196 if(l1>=0 && l1<Nz) Rdis(l1) += R(l);
197 }
198 }
199 // On remplit eventuellement les matrices 2D
200 if(do2d && M2d.Size()>0)
201 for(int l=0;l<Nz;l++) {M2d(j,l) = R(l); M2dc(j,l) = Rdis(l);}
202 // On remplit le cube avec le champ R redshift distordu
203 if(fluct3d) for(int l=0;l<Nz;l++) (*rgen)(l,j,i) += Rdis(l);
204 // Calcul eventuel de la fonction de correlation R*V
205 if(npixcor>0) {
206 for(long l1=0;l1<Nz;l1++) {
207 for(long l2=max(0L,l1-npixcor);l2<min(Nz,l1+npixcor);l2++) {
208 int lc = imil+(l2-l1);
209 Ksirr(lc) += R(l1)*R(l2);
210 Ksirrc(lc) += Rdis(l1)*R(l2);
211 Ksirv(lc) += R(l1)*V(l2);
212 Ksirvc(lc) += Rdis(l1)*V(l2);
213 nKsi(lc)++;
214 }
215 }
216 }
217 // Cross-power spectrum computation
218 if(dopk) {
219 fftserv.FFTForward(V,FV);
220 int nf = FV.Size();
221 if(pkvr.Size()<=0) {
222 cout<<"...Create vector for cross-power spectrum computation"<<endl;
223 pkvr.ReSize(nf); pkvr = complex<r_8>(0.);
224 pkr.ReSize(nf); pkr = 0.;
225 pkrc.ReSize(nf); pkrc = 0.;
226 }
227 fftserv.FFTForward(R,FR);
228 for(int l=0;l<nf;l++) {
229 pkvr(l) += FV(l)*conj(FR(l));
230 pkr(l) += norm(FR(l));
231 }
232 fftserv.FFTForward(Rdis,FRdis);
233 for(int l=0;l<nf;l++) pkrc(l) += norm(FRdis(l));
234 npk++;
235 }
236 nlosread++;
237 }
238 if(do2d && M2d.Size()>0) {
239 char str[64];
240 sprintf(str,"mx_%d",i);
241 pos.PutObject(M2d,str);
242 sprintf(str,"mxc_%d",i);
243 pos.PutObject(M2dc,str);
244 }
245 }
246
247 cout<<"Number of processed los: "<<nlosread<<" / "<<Nx*Ny<<endl;
248 dvlcor("nlosread") = (int_4)nlosread;
249 if(hmpc.NEntries()>0) {
250 hmpc.Show();
251 pos.PutObject(hmpc,"hmpc");
252 }
253 if(Ksirr.Size()>0) {
254 for(int l=0;l<Ksirr.Size();l++) if(nKsi(l)>0) {
255 Ksirr(l) /= nKsi(l);
256 Ksirrc(l) /= nKsi(l);
257 Ksirv(l) /= nKsi(l);
258 Ksirvc(l) /= nKsi(l);
259 }
260 pos.PutObject(Ksirr,"ksirr");
261 pos.PutObject(Ksirrc,"ksirrc");
262 pos.PutObject(Ksirv,"ksirv");
263 pos.PutObject(Ksirvc,"ksirvc");
264 pos.PutObject(nKsi,"nksi");
265 }
266 if(npk>0) {
267 pkvr /= (double)npk;
268 pkr /= (double)npk;
269 pkrc /= (double)npk;
270 pos.PutObject(pkvr,"pkvr");
271 pos.PutObject(pkr,"pkr");
272 pos.PutObject(pkrc,"pkrc");
273 }
274 PrtTim(">>>> End filling redshift distorted cube");
275
276 // --- Fourier transform 3D cube and compute 1D and 2D power spectra
277 if(fluct3d) {
278 cout<<">>> Fourier transform 3D cube and compute 1D and 2D power spectra"<<endl;
279 // do the FFT for spectrum analysis
280 fluct3d->ReComputeFourier();
281 PrtTim(">>>> End ReComputing spectrum");
282
283 // Compute 1D spectrum
284 cout<<endl<<"\n--- Computing final 1D spectrum"<<endl;
285 double dkmin = fluct3d->GetKincMin(), knyqmax = fluct3d->GetKmax();
286 long nherr = long(knyqmax/dkmin+0.5);
287 cout<<"\nFor HistoErr: d="<<dkmin<<" max="<<knyqmax<<" n="<<nherr<<endl;
288 HistoErr hpkrec(0.,knyqmax,nherr); hpkrec.Zero();
289 hpkrec.ReCenterBin(); hpkrec.Show();
290 fluct3d->ComputeSpectrum(hpkrec);
291 pos.PutObject(hpkrec,"hpkrec");
292 PrtTim(">>>> End Computing final spectrum");
293
294 // Compute 2D spectrum
295 cout<<"\n--- Computing final 2D spectrum"<<endl;
296 double dktmin = fluct3d->GetKTincMin(), ktnyqmax = fluct3d->GetKTmax();
297 double dkzmin = fluct3d->GetKinc()[2], kznyqmax = fluct3d->GetKnyq()[2];
298 long nherrt = long(ktnyqmax/dktmin+0.5), nherrz = long(kznyqmax/dkzmin+0.5);
299 cout<<"For Histo2DErr: d="<<dktmin<<","<<dkzmin
300 <<" max="<<ktnyqmax<<","<<kznyqmax<<" n="<<nherrt<<","<<nherrz<<endl;
301 Histo2DErr hpkrec2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
302 hpkrec2.ReCenterBin(); hpkrec2.Zero(); hpkrec2.Show();
303 fluct3d->ComputeSpectrum2D(hpkrec2);
304 pos.PutObject(hpkrec2,"hpkrec2");
305 PrtTim(">>>> End Computing final 2D spectrum");
306 }
307
308 // --- end of job, write objects in ppf
309 pos.PutObject(dvlcor,"dvlcor");
310 if(fluct3d) delete fluct3d;
311
312 //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
313 } catch (PException& exc) {
314 cerr<<"cmvrvloscor.cc catched PException"<<exc.Msg()<<endl;
315 return 77;
316 } catch (std::exception& sex) {
317 cerr << "cmvrvloscor.cc std::exception :"
318 << (string)typeid(sex).name() << "\n msg= "
319 << sex.what() << endl;
320 return 78;
321 } catch (...) {
322 cerr << "cmvrvloscor.cc catched unknown (...) exception " << endl;
323 return 79;
324 }
325 //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
326
327 return 0;
328}
329
330/*
331openppf cmvrvloscor.ppf
332
333disp hmpc
334
335# cross-correlation
336disp nksi
337set imil ${dvlcor.imil}
338n/plot ksirv.val%n-${imil} ! ! "nsta cpts"
339n/plot ksirvc.val%n-${imil} ! ! "nsta cpts same red"
340
341n/plot ksirr.val%n-${imil} ! ! "nsta cpts"
342n/plot ksirrc.val%n-${imil} ! ! "nsta cpts same red"
343
344# cross-power spectrum
345n/plot pkr.val%n ! ! "nsta cpts logx"
346n/plot pkrc.val%n ! ! "nsta cpts logx same red"
347
348n/plot pkvr.val%n ! ! "nsta cpts logx"
349
350# reconstructed 1D power spectrum
351n/plot hpkrec.val%x x>0 ! "nsta cpts logx"
352
353# reconstructed 2D power spectrum
354imag hpkrec2
355addoval 0 0 0.05 0.05 "green" false
356addoval 0 0 0.1 0.1 "green" false
357addoval 0 0 0.25 0.25 "green" false
358addoval 0 0 0.5 0.5 "green" false
359x = ${hpkrec2.xmax} / 2.
360addoval 0 0 $x $x "green" false
361x = ${hpkrec2.ymax} / 2.
362addoval 0 0 $x $x "green" false
363
364# proj selon kT (black), selon kZ (red)
365n/plot hpkrec2.val%sqrt(x*x+y*y) ! ! "nsta crossmarker3 logx"
366
367defscript tom
368del m2
369c++exec TMatrix<r_8> m2(hpkrec2.NBinX(),hpkrec2.NBinY()); \
370for(int i=0;i<hpkrec2.NBinX();i++) \
371for(int j=0;j<hpkrec2.NBinY();j++) m2(i,j) = hpkrec2(i,j); \
372KeepObj(m2); cout<<"OK"<<endl;
373endscript
374
375# les matrices 2D
376set n 0
377disp mx_$n
378newwin
379disp mxc_$n
380 */
Note: See TracBrowser for help on using the repository browser.