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

#####################################################################################
####  Commands to run the different programs to produce foreground maps 
####  and compute radio-source subtracted P(k)
#####################################################################################

### Cube definition in file cubedef.h 

### Step 1/ Produce an LSS data cube  with appropriate size and redshift using SimLSS
# 1.a/ Run SimLSS
csh> ~/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
# 1.b/ To run SimLSS with GSM map parameters (DeltaFreq=500 MHz)
csh> ~/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

# 1.c/ To run SimLSS with GSM map parametersand 40x40 deg maps (DeltaFreq=500 MHz)
csh> ~/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

# 1.c/ Change the X and Z axis of the cube to adapt it to RadioBeam package convention 
#  SimLSS output : the radial (redshift) direction along X axis of the cube (TArray) 
#  RadioBeam cubes : the radial (redshift) direction along Z axis of the cube (TArray) 
#  Execucte the following script in spiapp :

csh> cat > racube.pic
set f lssz060
readfits ${f}_r.fits
rename ${f}_r map
print map 
c++exec \
  TArray<r_4> omap(map.SizeY(),map.SizeZ(),map.SizeX()-2 ); \
  for(sa_size_t i=0;i<omap.SizeX();i++) \
    for(sa_size_t j=0;j<omap.SizeY();j++) \
      for(sa_size_t k=0;k<omap.SizeZ();k++) \
        omap(i,j,k)=map(k+1,i,j);  \
  KeepObj(omap); 

rename omap lsscube
print lsscube
# expmeansig lsscube val
saveppf lsscube lsscubez060.ppf 

csh> spiapp -term -exec racube.pic 

#### Cube LSS 40x30 deg (3')  @ z=0.6 ( lsscubez060.ppf ) 
#### -> Size= 122880000 Mean=-7.01664e-05 Sigma=2.53016 Min=-13.7439 Max=14.4648

## Step 2/ Produce synchrotron and radio source sky cubes  (cube unit is Temparature- Kelvin) 
# 2.a/ Synchrotron map from HASLAM 400 MHz map
csh> ./Objs/syncube syncmap_eq.fits syncube.ppf syncmap.ppf
# 2.b/ radio source cube from NVSS catalog
csh> ./Objs/srcat2cube -nvss nvss.fits nvsscube.ppf nvssmap.ppf
# Or from the north20 catalog :
csh> ./Objs/srcat2cube -north20 north20cm.fits north20cube.ppf north20map.ppf

# 2.c/ Add the two cubes using the following spiapp script 
csh> cat > sumcubes.pic
openppf syncube.ppf
openppf radsrccube.ppf
# expmeansig syncube val 
# expmeansig nvsscube val
c++exec TArray<r_4> fgndcube = syncube+radsrccube; KeepObj(fgndcube);
print fgndcube
# expmeansig fgndcube val
saveppf fgndcube fgndcube.ppf

csh> spiapp -term -exec sumcubes.pic 

#### syncube:Mean= 1.8101  Sigma= 0.326538  Min= 0.857019 Max= 3.58987
#### nvsscube: Mean= 1.95073  Sigma= 1.68515  Min= 0.857019 Max= 428.398
#### fgndcube=syncube+nvsscube: Mean= 0.140623  Sigma= 1.65068  Min= 0 Max= 426.559
#### north20: fgndcube_north:

## Step 2.b/ Produce foreground cube from GSM 
csh> ./Objs/gsm2cube ../Catalogs/GSM/ 1 256 fgndcube_gsm.ppf

## Step 3/ Apply lobe (50 meter diameter array) effect on foreground cube and LSS cube
csh> set ddish=55.
csh> set ddishcor=55.
csh> ./Objs/applobe $ddish fgndcube.ppf fgndcube_lobe.ppf
csh> ./Objs/applobe -fib $ddish fgndcube.ppf fgndcube_flobe.ppf
csh> ./Objs/applobe $ddish lsscube.ppf lsscube_lobe.ppf
csh> ./Objs/applobe -fib $ddish lsscube.ppf lsscube_flobe.ppf
## 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)  
csh> ./Objs/applobe $ddish lsscube_lobe.ppf lsscube_corlobe.ppf $ddishcor
csh> ./Objs/applobe $ddish fgndcube_lobe.ppf fgndcube_corlobe.ppf $ddishcor 

## Step 3.c/ Apply lobe (Filled 11x11 5m dishes array) effect on foreground cube and LSS cube
csh> ./Objs/applobe repf11x11.ppf fgndcube.ppf fgndcube_lobe.ppf
csh> ./Objs/applobe repf11x11.ppf lsscube.ppf lsscube_lobe.ppf
## 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) 
csh> ./Objs/applobe repf11x11.ppf lsscube_lobe.ppf lsscube_corlobe.ppf $ddishcor
csh> ./Objs/applobe repf11x11.ppf fgndcube_lobe.ppf fgndcube_corlobe.ppf $ddishcor

### Step 4/ Compute power spectra 
## 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)
## Foreground maps are in temperature 
## Noise fluctuations Sigma^2 ~ T_sys^2 / t_obs * DeltaFreq
## Tsys ~ 50 K , DeltaFreq ~ 0.5 MHz , t_obs ~ 1 day ~ 80 000 s.
## sigma_noise ~ 0.25 mK -> 3 mK  
# 4.a/ LSS power spectrum  without noise
csh> ./Objs/calcpk lsscube.ppf lsspk.ppf 0.13
# and with noise  
csh> ./Objs/calcpk lsscube.ppf lsspkwn.ppf 0.13  3
#  with the lobe effect 
csh> ./Objs/calcpk lsscube_lobe.ppf lsspklobe.ppf 0.13
csh> ./Objs/calcpk lsscube_flobe.ppf lsspkflobe.ppf 0.13
csh> ./Objs/calcpk lsscube_lobe.ppf lsspklobewn.ppf 0.13 3  
csh> ./Objs/calcpk lsscube_corlobe.ppf lsspkcorlobe.ppf 0.13  

# 4.b/ Foreground power spectrum 
csh> ./Objs/calcpk fgndcube.ppf fgndpk.ppf 1000
csh> ./Objs/calcpk fgndcube_lobe.ppf fgndpklobe.ppf 1000
csh> ./Objs/calcpk fgndcube_flobe.ppf fgndpkflobe.ppf 1000
csh> ./Objs/calcpk fgndcube_corlobe.ppf fgndpkcorlobe.ppf 1000

# 4.c/ Extract LSS P(k) from Foreground+LSS+noise , after cleaning/subtraction without beam 
csh> set beamdesc=repf11x11.ppf
csh> set ddishcor=55.
csh> set noiselev=1.
csh> ./Objs/calcpk2 lsscube.ppf 0.13 fgndcube.ppf 1000 subpk.ppf $noiselev $beamdesc 0. 0. P2
# 4.d / Extract LSS P(k) from Foreground+LSS+noise and beam effect, without beam correction  
csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.13 fgndcube_lobe.ppf 1000 subpklobe.ppf $noiselev $beamdesc 0. 0. P2 
# 4.e / Extract LSS P(k) from Foreground+LSS+noise and beam effect - correcting to a beam of Diam= $ddishcor
csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.13 fgndcube_lobe.ppf 1000 subpkcorlobe.ppf $noiselev $beamdesc $ddishcor 0. P2 reclsscorlobe.ppf 
#  Or using a linear fit for foreground subtraction (old version)
csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.13 fgndcube_lobe.ppf 1000 subpkcorlobep1.ppf $noiselev $beamdesc $ddishcor 0. P2 
# 4.f / Estimate residual noise from Foreground removal : 
csh> ./Objs/calcpk2 lsscube.ppf 0. fgndcube_lobe.ppf 1000 residcorlobe.ppf  $noiselev $beamdesc $ddishcor 0. P2
csh> ./Objs/calcpk2 lsscube.ppf 0. fgndcube_lobe.ppf 1000 residnocor.ppf $noiselev $beamdesc 0. 0. P2

### Step 5 / Check the results using spiapp 
setaxesatt 'font=helvetica,bold,16 fixedfontsize minorticks'
delobjs *
openppf fgndpk.ppf 
openppf fgndpklobe.ppf 
openppf fgndpkflobe.ppf 
openppf fgndpkcorlobe.ppf 
openppf lsspk.ppf 
openppf lsspklobe.ppf 
openppf lsspkcorlobe.ppf 

openppf lsspkflobe.ppf 
openppf lsspklobewn.ppf 
openppf subpklobe.ppf 

openppf subpkcorlobe.ppf 

openppf residcorlobe.ppf 
openppf residnocor.ppf
# openppf subpknolssnocor.ppf

disp lsspk 'logx logy nsta xylimits=0.01,2.,4e-11,8e-6 gold'
disp lsspklobe 'same nsta orange'
disp lsspklobewn 'same nsta siennared'
settitle ' Pk[LSS] - without normalisation' ' ' 'font=helvetica,bold,16 black' 


disp fgndpk 'logx logy nsta xylimits=0.01,2.,1e-10,1. navyblue'
disp fgndpklobe 'same nsta blue'
disp fgndpkcorlobe 'same nsta skyblue'
disp lsspk 'same nsta gold'
disp lsspkflobe 'same nsta yellow'
disp subpkcorlobe 'same nsta red' 
disp residcorlobe 'same nsta green'
disp residnocor 'same nsta forestgreen'

disp lsspklobewn 'same nsta siennared' 
# settitle 'Pk[LSS] , Pk[Foreground] and lobe effect (Dish D=50 m)' ' ' 'font=helvetica,bold,18'
settitle 'Pk[LSS] , Pk[Foreground=GSM] and lobe effect (Dish D=50 m)' ' ' 'font=helvetica,bold,18'

set lines ( 'Pk[Foreground]'  'Pk[fgnd]*Lobe' 'Pk[fgnd]*Lobe/Corrected' 'Pk[LSS]' 'Pk[LSS]*Lobe+Noise' )
set cols ( navyblue blue skyblue gold siennared ) 
textdrawer lines cols 'font=helvetica,bold,16 frame'


disp lsspk 'logx logy nsta xylimits=0.005,2.,4e-9,4e-5 gold'
disp lsspklobewn 'same nsta siennared' 
disp subpkcorlobe 'same nsta red'
disp subpknolss 'same nsta green'

#  Calcul du volume total en Mpc^3
set VOL (1.9*1.9*2.8*800*600*256)
# set VOL (1.9*1.9*2.8*1800*600*256)
plot2d 
# plot2d lsspk x val*$VOL 1 'logx logy nsta xylimits=0.01,2.,10.,1e4 cpts marker=box,5 gold'
plot2d lsspklobewn x val*$VOL 1 'same nsta cpts marker=box,5 siennared' 
plot2d subpkcorlobe x val*$VOL 1 'same nsta cpts marker=box,5 red' 
plot2d subpklobe x val*$VOL 1 'same nsta cpts marker=box,5 blueviolet' 
plot2d subpknolss x val*$VOL 1 'same nsta cpts marker=box,5 green' 

# settitle 'Recovered Pk[LSS] In=LSS+(GSM) (D=50 m)' ' ' 'font=helvetica,bold,18'
settitle 'Recovered Pk[LSS] In=LSS+(Haslam+North20cm) (D=50 m)' ' ' 'font=helvetica,bold,18'
setaxelabels 'k (Mpc^-1)   h=0.7' 'P(k)    (mK^2 Mpc^3)' 'font=helvetica,bolditalic,16' 
set lines ( 'Pk[LSS]'  'Pk[LSS*lobe+noise]' 'Pk[ExtractedLSS]' 'Pk[ExtLSS,NoBeamCor]' 'Pk[residual,NoLSS]' )
set cols ( gold siennared red blueviolet green ) 
textdrawer lines cols 'font=helvetica,bold,16 frame'

plot2d fgndpk x val*$VOL 1 'logx logy  xylimits=0.01,1.,1.,1e10 nsta cpts marker=box,5 black'
plot2d fgndpkflobe x val*$VOL 1 ' nsta cpts marker=circle,5 navyblue same'
plot2d fgndpklobe x val*$VOL 1 ' nsta cpts marker=circle,5 blue same'

plot2d lsspk x val*$VOL 1 ' nsta cpts marker=box,5 red same'
plot2d lsspkflobe x val*$VOL 1 ' nsta cpts marker=circle,5 orange same'
plot2d lsspklobe x val*$VOL 1 ' nsta cpts marker=circle,5 yellow same'

c++exec \
  Histo lsspkratioA = subpkcorlobe/lsspk; KeepObj(lsspkratioA); \
  Histo lsspkratioB = subpkcorlobe/lsspkflobe; KeepObj(lsspkratioB); 

plot2d lsspk x val*$VOL 1 'logx logy nsta xylimits=0.01,2.,10.,1e4 cpts marker=box,5 gold'
plot2d lsspklobewn x val*$VOL 1 'same nsta cpts marker=box,5 red' 
plot2d subpknolss x val*$VOL 1 'same nsta cpts marker=box,5 green' 
plot2d subpknolssnocor x val*$VOL 1 'same nsta cpts marker=box,5 magenta' 
setaxelabels 'k (Mpc^-1)   h=0.7' 'P(k)    (mK^2 Mpc^3)' 'font=helvetica,bolditalic,16' 
settitle 'Recovered Pk[LSS] and residual systematics' ' ' 'font=helvetica,bold,18'
set lines ( 'Pk[LSS]'  'Pk[LSS*lobe+noise]' 'Pk[residual,NoLSS]' 'Pk[residual,NoLSS,NoBeamCorrection]' )
set cols ( gold red green magenta ) 
textdrawer lines cols 'font=helvetica,bold,16 frame'