source: Sophya/trunk/Cosmo/RadioBeam/subtractradsrc.cmd@ 4045

Last change on this file since 4045 was 3986, checked in by ansari, 14 years ago

modification rapport maxi a appliquer lors des corrections de beams, Reza 05/05/2011

File size: 10.4 KB
Line 
1#####################################################################################
2#### Commands to run the different programs to produce foreground maps
3#### and compute radio-source subtracted P(k)
4#####################################################################################
5
6### Cube definition in file cubedef.h
7
8### Step 1/ Produce an LSS data cube with appropriate size and redshift using SimLSS
9# 1.a/ Run SimLSS
10csh> ~/Objs/exe/cmvginit3df -a -1 -2 -C -G 0. -F 0 -x 360,3 -y 360,3 -z 256,1.5 -Z 0.56 -8 1. -n 10000 -O 0,2 -o lssz056 -T 2
11# 1.b/ To run SimLSS with GSM map parameters (DeltaFreq=500 MHz)
12csh> ~/Objs/exe/cmvginit3df -a -1 -2 -C -G 0. -F 0 -x 360,3 -y 360,3 -z 256,3 -Z 0.60 -8 1. -n 10000 -O 0,2 -o lssz060 -T 2
13
14# 1.c/ To run SimLSS with GSM map parametersand 40 (phi/alpha) x 30 (theta/delta) deg maps (DeltaFreq=500 MHz) @ z=0.6
15csh> ~/Objs/exe/cmvginit3df -a -1 -2 -C -G 0. -F 0 -x 600,1.9 -y 800,1.9 -z 256,2.8 -Z 0.60 -8 1. -n 10000 -O 0,2 -o lssz060 -T 2
16
17# 1.d/ To run SimLSS with GSM map parametersand 90x30 deg maps (DeltaFreq=500 MHz) @ z=1 [ 90deg-> phi/alpha, 30deg -> theta/delta]
18csh> ~/Objs/exe/cmvginit3df -a -1 -2 -C -G 0. -F 0 -x 600,2.9 -y 1800,2.9 -z 256,3.5 -Z 1.0 -8 1. -n 10000 -O 0,2 -o lssz100 -T 2
19
20csh> ~/Objs/exe/cmvginit3df -a -1 -2 -C -G 0. -F 0 -x 600,3.8 -y 1800,3.8 -z 256,4.2 -Z 1.5 -8 1. -n 10000 -O 0,2 -o lssz150 -T 2
21
22# 1.c/ Change the X and Z axis of the cube to adapt it to RadioBeam package convention
23# SimLSS output : the radial (redshift) direction along X axis of the cube (TArray)
24# RadioBeam cubes : the radial (redshift) direction along Z axis of the cube (TArray)
25# Execucte the following script in spiapp :
26
27csh> cat > racube.pic
28set f lssz060
29readfits ${f}_r.fits
30rename ${f}_r map
31print map
32c++exec \
33 TArray<r_4> omap(map.SizeY(),map.SizeZ(),map.SizeX()-2 ); \
34 for(sa_size_t i=0;i<omap.SizeX();i++) \
35 for(sa_size_t j=0;j<omap.SizeY();j++) \
36 for(sa_size_t k=0;k<omap.SizeZ();k++) \
37 omap(i,j,k)=map(k+1,i,j); \
38 KeepObj(omap);
39
40rename omap lsscube
41print lsscube
42# expmeansig lsscube val
43saveppf lsscube lsscubez060.ppf
44
45csh> spiapp -term -exec racube.pic
46
47#### Cube LSS 40x30 deg (3') @ z=0.6 ( lsscubez060.ppf )
48#### -> Size= 122880000 Mean=-7.01664e-05 Sigma=2.53016 Min=-13.7439 Max=14.4648
49
50## Step 2/ Produce synchrotron and radio source sky cubes (cube unit is Temparature- Kelvin)
51# 2.a/ Synchrotron map from HASLAM 400 MHz map
52csh> ./Objs/syncube syncmap_eq.fits syncube.ppf syncmap.ppf
53# 2.b/ radio source cube from NVSS catalog
54csh> ./Objs/srcat2cube -nvss nvss.fits nvsscube.ppf nvssmap.ppf
55# Or from the north20 catalog :
56csh> ./Objs/srcat2cube -north20 north20cm.fits north20cube.ppf north20map.ppf
57
58# 2.c/ Add the two cubes using the following spiapp script
59csh> cat > sumcubes.pic
60openppf syncube.ppf
61openppf radsrccube.ppf
62# expmeansig syncube val
63# expmeansig nvsscube val
64c++exec TArray<r_4> fgndcube = syncube+radsrccube; KeepObj(fgndcube);
65print fgndcube
66# expmeansig fgndcube val
67saveppf fgndcube fgndcube.ppf
68
69csh> spiapp -term -exec sumcubes.pic
70
71#### syncube:Mean= 1.8101 Sigma= 0.326538 Min= 0.857019 Max= 3.58987
72#### nvsscube: Mean= 1.95073 Sigma= 1.68515 Min= 0.857019 Max= 428.398
73#### fgndcube=syncube+nvsscube: Mean= 0.140623 Sigma= 1.65068 Min= 0 Max= 426.559
74#### north20: fgndcube_north:
75
76## Step 2.b/ Produce foreground cube from GSM
77csh> ./Objs/gsm2cube ../Catalogs/GSM/ 1 256 fgndcube_gsm.ppf
78
79## Step 3/ Apply lobe (50 meter diameter array) effect on foreground cube and LSS cube
80csh> set ddish=55.
81csh> set ddishcor=55.
82csh> ./Objs/applobe $ddish fgndcube.ppf fgndcube_lobe.ppf
83csh> ./Objs/applobe -fib $ddish fgndcube.ppf fgndcube_flobe.ppf
84csh> ./Objs/applobe $ddish lsscube.ppf lsscube_lobe.ppf
85csh> ./Objs/applobe -fib $ddish lsscube.ppf lsscube_flobe.ppf
86## Step 3.b/ Correct for the lobe effect by bringing all to the beam of Diam/Lambda = 150 (55 m @ z=0.7 - 820 MHz)
87csh> ./Objs/applobe $ddish lsscube_lobe.ppf lsscube_corlobe.ppf $ddishcor
88csh> ./Objs/applobe $ddish fgndcube_lobe.ppf fgndcube_corlobe.ppf $ddishcor
89
90## Step 3.c/ Apply lobe (Filled 11x11 5m dishes array) effect on foreground cube and LSS cube
91csh> ./Objs/applobe repf11x11.ppf fgndcube.ppf fgndcube_lobe.ppf
92csh> ./Objs/applobe repf11x11.ppf lsscube.ppf lsscube_lobe.ppf
93## Step 3.d/ Correct for the lobe effect by bringing all to the beam of Diam/Lambda = 150 (55 m @ z=0.7 - 820 MHz)
94csh> ./Objs/applobe repf11x11.ppf lsscube_lobe.ppf lsscube_corlobe.ppf $ddishcor
95csh> ./Objs/applobe repf11x11.ppf fgndcube_lobe.ppf fgndcube_corlobe.ppf $ddishcor
96
97### Step 4/ Compute power spectra
98## mass to temperature converion factor CT21 ~= 0.21 mK for gHI=2% , 0.11 for gHI=1% , 0.13 for gHI=0.008x(1+0.6)
99## Foreground maps are in temperature
100## Noise fluctuations Sigma^2 ~ T_sys^2 / t_obs * DeltaFreq
101## Tsys ~ 50 K , DeltaFreq ~ 0.5 MHz , t_obs ~ 1 day ~ 80 000 s.
102## sigma_noise ~ 0.25 mK -> 3 mK
103# 4.a/ LSS power spectrum without noise
104csh> ./Objs/calcpk lsscube.ppf lsspk.ppf 0.13
105# and with noise
106csh> ./Objs/calcpk lsscube.ppf lsspkwn.ppf 0.13 3
107# with the lobe effect
108csh> ./Objs/calcpk lsscube_lobe.ppf lsspklobe.ppf 0.13
109csh> ./Objs/calcpk lsscube_flobe.ppf lsspkflobe.ppf 0.13
110csh> ./Objs/calcpk lsscube_lobe.ppf lsspklobewn.ppf 0.13 3
111csh> ./Objs/calcpk lsscube_corlobe.ppf lsspkcorlobe.ppf 0.13
112
113# 4.b/ Foreground power spectrum
114csh> ./Objs/calcpk fgndcube.ppf fgndpk.ppf 1000
115csh> ./Objs/calcpk fgndcube_lobe.ppf fgndpklobe.ppf 1000
116csh> ./Objs/calcpk fgndcube_flobe.ppf fgndpkflobe.ppf 1000
117csh> ./Objs/calcpk fgndcube_corlobe.ppf fgndpkcorlobe.ppf 1000
118
119# 4.c/ Extract LSS P(k) from Foreground+LSS+noise , after cleaning/subtraction without beam
120csh> set beamdesc=repf11x11.ppf
121csh> set ddishcor=55.
122csh> set noiselev=1.
123csh> ./Objs/calcpk2 lsscube.ppf 0.13 fgndcube.ppf 1000 subpk.ppf $noiselev $beamdesc 0. 0. P2
124# 4.d / Extract LSS P(k) from Foreground+LSS+noise and beam effect, without beam correction
125csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.13 fgndcube_lobe.ppf 1000 subpklobe.ppf $noiselev $beamdesc 0. 0. P2
126# 4.e / Extract LSS P(k) from Foreground+LSS+noise and beam effect - correcting to a beam of Diam= $ddishcor
127csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.13 fgndcube_lobe.ppf 1000 subpkcorlobe.ppf $noiselev $beamdesc $ddishcor 0. P2 reclsscorlobe.ppf
128# Or using a linear fit for foreground subtraction (old version)
129csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.13 fgndcube_lobe.ppf 1000 subpkcorlobep1.ppf $noiselev $beamdesc $ddishcor 0. P2
130# 4.f / Estimate residual noise from Foreground removal :
131csh> ./Objs/calcpk2 lsscube.ppf 0. fgndcube_lobe.ppf 1000 residcorlobe.ppf $noiselev $beamdesc $ddishcor 0. P2
132csh> ./Objs/calcpk2 lsscube.ppf 0. fgndcube_lobe.ppf 1000 residnocor.ppf $noiselev $beamdesc 0. 0. P2
133
134### Step 5 / Check the results using spiapp
135setaxesatt 'font=helvetica,bold,16 fixedfontsize minorticks'
136delobjs *
137openppf fgndpk.ppf
138openppf fgndpklobe.ppf
139openppf fgndpkflobe.ppf
140openppf fgndpkcorlobe.ppf
141openppf lsspk.ppf
142openppf lsspklobe.ppf
143openppf lsspkcorlobe.ppf
144
145openppf lsspkflobe.ppf
146openppf lsspklobewn.ppf
147openppf subpklobe.ppf
148
149openppf subpkcorlobe.ppf
150
151openppf residcorlobe.ppf
152openppf residnocor.ppf
153# openppf subpknolssnocor.ppf
154
155disp lsspk 'logx logy nsta xylimits=0.01,2.,4e-11,8e-6 gold'
156disp lsspklobe 'same nsta orange'
157disp lsspklobewn 'same nsta siennared'
158settitle ' Pk[LSS] - without normalisation' ' ' 'font=helvetica,bold,16 black'
159
160
161disp fgndpk 'logx logy nsta xylimits=0.01,2.,1e-10,1. navyblue'
162disp fgndpklobe 'same nsta blue'
163disp fgndpkcorlobe 'same nsta skyblue'
164disp lsspk 'same nsta gold'
165disp lsspkflobe 'same nsta yellow'
166disp subpkcorlobe 'same nsta red'
167disp residcorlobe 'same nsta green'
168disp residnocor 'same nsta forestgreen'
169
170disp lsspklobewn 'same nsta siennared'
171# settitle 'Pk[LSS] , Pk[Foreground] and lobe effect (Dish D=50 m)' ' ' 'font=helvetica,bold,18'
172settitle 'Pk[LSS] , Pk[Foreground=GSM] and lobe effect (Dish D=50 m)' ' ' 'font=helvetica,bold,18'
173
174set lines ( 'Pk[Foreground]' 'Pk[fgnd]*Lobe' 'Pk[fgnd]*Lobe/Corrected' 'Pk[LSS]' 'Pk[LSS]*Lobe+Noise' )
175set cols ( navyblue blue skyblue gold siennared )
176textdrawer lines cols 'font=helvetica,bold,16 frame'
177
178
179disp lsspk 'logx logy nsta xylimits=0.005,2.,4e-9,4e-5 gold'
180disp lsspklobewn 'same nsta siennared'
181disp subpkcorlobe 'same nsta red'
182disp subpknolss 'same nsta green'
183
184# Calcul du volume total en Mpc^3
185set VOL (1.9*1.9*2.8*800*600*256)
186# set VOL (1.9*1.9*2.8*1800*600*256)
187plot2d
188# plot2d lsspk x val*$VOL 1 'logx logy nsta xylimits=0.01,2.,10.,1e4 cpts marker=box,5 gold'
189plot2d lsspklobewn x val*$VOL 1 'same nsta cpts marker=box,5 siennared'
190plot2d subpkcorlobe x val*$VOL 1 'same nsta cpts marker=box,5 red'
191plot2d subpklobe x val*$VOL 1 'same nsta cpts marker=box,5 blueviolet'
192plot2d subpknolss x val*$VOL 1 'same nsta cpts marker=box,5 green'
193
194# settitle 'Recovered Pk[LSS] In=LSS+(GSM) (D=50 m)' ' ' 'font=helvetica,bold,18'
195settitle 'Recovered Pk[LSS] In=LSS+(Haslam+North20cm) (D=50 m)' ' ' 'font=helvetica,bold,18'
196setaxelabels 'k (Mpc^-1) h=0.7' 'P(k) (mK^2 Mpc^3)' 'font=helvetica,bolditalic,16'
197set lines ( 'Pk[LSS]' 'Pk[LSS*lobe+noise]' 'Pk[ExtractedLSS]' 'Pk[ExtLSS,NoBeamCor]' 'Pk[residual,NoLSS]' )
198set cols ( gold siennared red blueviolet green )
199textdrawer lines cols 'font=helvetica,bold,16 frame'
200
201plot2d fgndpk x val*$VOL 1 'logx logy xylimits=0.01,1.,1.,1e10 nsta cpts marker=box,5 black'
202plot2d fgndpkflobe x val*$VOL 1 ' nsta cpts marker=circle,5 navyblue same'
203plot2d fgndpklobe x val*$VOL 1 ' nsta cpts marker=circle,5 blue same'
204
205plot2d lsspk x val*$VOL 1 ' nsta cpts marker=box,5 red same'
206plot2d lsspkflobe x val*$VOL 1 ' nsta cpts marker=circle,5 orange same'
207plot2d lsspklobe x val*$VOL 1 ' nsta cpts marker=circle,5 yellow same'
208
209c++exec \
210 Histo lsspkratioA = subpkcorlobe/lsspk; KeepObj(lsspkratioA); \
211 Histo lsspkratioB = subpkcorlobe/lsspkflobe; KeepObj(lsspkratioB);
212
213plot2d lsspk x val*$VOL 1 'logx logy nsta xylimits=0.01,2.,10.,1e4 cpts marker=box,5 gold'
214plot2d lsspklobewn x val*$VOL 1 'same nsta cpts marker=box,5 red'
215plot2d subpknolss x val*$VOL 1 'same nsta cpts marker=box,5 green'
216plot2d subpknolssnocor x val*$VOL 1 'same nsta cpts marker=box,5 magenta'
217setaxelabels 'k (Mpc^-1) h=0.7' 'P(k) (mK^2 Mpc^3)' 'font=helvetica,bolditalic,16'
218settitle 'Recovered Pk[LSS] and residual systematics' ' ' 'font=helvetica,bold,18'
219set lines ( 'Pk[LSS]' 'Pk[LSS*lobe+noise]' 'Pk[residual,NoLSS]' 'Pk[residual,NoLSS,NoBeamCorrection]' )
220set cols ( gold red green magenta )
221textdrawer lines cols 'font=helvetica,bold,16 frame'
Note: See TracBrowser for help on using the repository browser.