1 | ############################################################## |
---|
2 | # Check chromaticity of DAB from calib spectra (1204-window) # |
---|
3 | ############################################################## |
---|
4 | |
---|
5 | set user "AST" |
---|
6 | set toppath "/sps/baoradio/AmasNancay/${user}" |
---|
7 | set source "Abell1205" |
---|
8 | set srclower "abell1205" |
---|
9 | set date "20110415" |
---|
10 | set cycle "8" |
---|
11 | set mode "Off" |
---|
12 | set path "${toppath}/${source}/${date}${srclower}/${mode}/calibcycle${cycle}" |
---|
13 | |
---|
14 | # |
---|
15 | # Power vs. run (calib monitoring) |
---|
16 | #----------------------------------- |
---|
17 | #shell ls ${path}/*.fits |
---|
18 | shell ls ${path}/*.fits | wc -l > pp |
---|
19 | vecfrascii mm pp |
---|
20 | vec2var mm nfiles |
---|
21 | |
---|
22 | set mean0 " " |
---|
23 | set mean1 " " |
---|
24 | for i 0:$nfiles |
---|
25 | openfits ${path}/medfiltmtx${i}.fits |
---|
26 | del g0 g1 s0 s1 n0 n1 |
---|
27 | objaoper /home/medfiltmtx${i} row 0 g0 |
---|
28 | objaoper /home/medfiltmtx${i} row 1 g1 |
---|
29 | c++exec sa_size_t n0=g0.NElts()-1; TVector<r_4>s0=g0(Range(1,n0)); KeepObj(s0); \ |
---|
30 | sa_size_t n1=g1.NElts()-1; TVector<r_4>s1=g1(Range(1,n1)); KeepObj(s1); |
---|
31 | m0 = ${s0.sum}/${s0.size} |
---|
32 | m1 = ${s1.sum}/${s1.size} |
---|
33 | set mean0 "${mean0} ${m0}" |
---|
34 | set mean1 "${mean1} ${m1}" |
---|
35 | end |
---|
36 | echo ${mean0} |
---|
37 | echo ${mean1} |
---|
38 | line2vec vmean0 ${mean0} |
---|
39 | line2vec vmean1 ${mean1} |
---|
40 | |
---|
41 | maxplot = ${vmean1.max}+5 |
---|
42 | newwin 1 1 |
---|
43 | graphicatt "xylimits=0,116,0,3" |
---|
44 | plot2d vmean1 n val 1 "red cpts nsta notit" |
---|
45 | plot2d vmean0 n val 1 "same blue cpts nsta notit" |
---|
46 | settitle "Mean calib $source ${date} ${mode} Ch 0 (blue) Ch 1 (red)" |
---|
47 | setaxelabels "File #" "I (a.u)" |
---|
48 | |
---|
49 | # |
---|
50 | # Identify baseline, DAB_1, DAB_2 from monitoring plot |
---|
51 | #------------------------------------------------------ |
---|
52 | ## Baseline |
---|
53 | del ifirst count nspec sum0 sum1 |
---|
54 | ifirst = 1 |
---|
55 | count = 0 |
---|
56 | for i 0:34 |
---|
57 | # echo "^^^^^^^^" $i $count |
---|
58 | count = $count+1 |
---|
59 | del g0 g1 |
---|
60 | objaoper /home/medfiltmtx${i} row 0 g0 |
---|
61 | objaoper /home/medfiltmtx${i} row 1 g1 |
---|
62 | if ( $ifirst == 1 ) then |
---|
63 | ifirst = 0 |
---|
64 | c++exec TVector<r_4> sum0=g0; KeepObj(sum0); \ |
---|
65 | TVector<r_4> sum1=g1; KeepObj(sum1); |
---|
66 | else |
---|
67 | c++exec sum0+=g0; sum1+=g1; |
---|
68 | endif |
---|
69 | end |
---|
70 | for i 70:$nfiles |
---|
71 | count = $count+1 |
---|
72 | # echo "^^^^^^^^" $i $count |
---|
73 | del g0 g1 |
---|
74 | objaoper /home/medfiltmtx${i} row 0 g0 |
---|
75 | objaoper /home/medfiltmtx${i} row 1 g1 |
---|
76 | c++exec sum0+=g0; sum1+=g1; |
---|
77 | end |
---|
78 | echo $count |
---|
79 | line2vec nspec $count |
---|
80 | del base0 base1 |
---|
81 | c++exec TVector<r_4> base0=sum0; KeepObj(base0); base0/=nspec[0]; \ |
---|
82 | TVector<r_4> base1=sum1; KeepObj(base1); base1/=nspec[0]; |
---|
83 | |
---|
84 | ## DAB Level 1 |
---|
85 | idablev1 = 36 |
---|
86 | idablev2 = 51 |
---|
87 | |
---|
88 | del ifirst count nspec sum0 sum1 |
---|
89 | ifirst = 1 |
---|
90 | count = 0 |
---|
91 | for i $idablev1:$idablev2 |
---|
92 | count = $count+1 |
---|
93 | echo "^^^^^^^^" $i $count |
---|
94 | del g0 g1 |
---|
95 | objaoper /home/medfiltmtx${i} row 0 g0 |
---|
96 | objaoper /home/medfiltmtx${i} row 1 g1 |
---|
97 | if ( $ifirst == 1 ) then |
---|
98 | ifirst = 0 |
---|
99 | c++exec TVector<r_4> sum0=g0; KeepObj(sum0); \ |
---|
100 | TVector<r_4> sum1=g1; KeepObj(sum1); |
---|
101 | else |
---|
102 | c++exec sum0+=g0; sum1+=g1; |
---|
103 | endif |
---|
104 | end |
---|
105 | echo $count |
---|
106 | line2vec nspec $count |
---|
107 | del dab01 dab11 |
---|
108 | c++exec TVector<r_4> dab01=sum0; KeepObj(dab01); dab01/=nspec[0]; \ |
---|
109 | TVector<r_4> dab11=sum1; KeepObj(dab11); dab11/=nspec[0]; |
---|
110 | |
---|
111 | ## DAB Level 2 |
---|
112 | idablev3 = 53 |
---|
113 | idablev4 = 67 |
---|
114 | |
---|
115 | del ifirst count nspec sum0 sum1 |
---|
116 | ifirst = 1 |
---|
117 | count = 0 |
---|
118 | for i $idablev3:$idablev4 |
---|
119 | count = $count+1 |
---|
120 | echo "^^^^^^^^" $i $count |
---|
121 | del g0 g1 |
---|
122 | objaoper /home/medfiltmtx${i} row 0 g0 |
---|
123 | objaoper /home/medfiltmtx${i} row 1 g1 |
---|
124 | if ( $ifirst == 1 ) then |
---|
125 | ifirst = 0 |
---|
126 | c++exec TVector<r_4> sum0 = g0; KeepObj(sum0); \ |
---|
127 | TVector<r_4> sum1 = g1; KeepObj(sum1); |
---|
128 | else |
---|
129 | c++exec sum0+=g0; sum1+=g1; |
---|
130 | endif |
---|
131 | end |
---|
132 | echo $count |
---|
133 | line2vec nspec $count |
---|
134 | del dab02 dab12 |
---|
135 | c++exec TVector<r_4> dab02=sum0; KeepObj(dab02); dab02/=nspec[0]; \ |
---|
136 | TVector<r_4> dab12=sum1; KeepObj(dab12); dab12/=nspec[0]; |
---|
137 | |
---|
138 | ## Plot results |
---|
139 | newwin 1 1 |
---|
140 | plot2d base0 (n/8192)*250+1250 val 1 "blue cpts notit nstat" |
---|
141 | plot2d base1 (n/8192)*250+1250 val 1 "same red cpts notit nstat" |
---|
142 | settitle "Mean baseline $source $date Ch0 blue Ch1 red" |
---|
143 | |
---|
144 | newwin 1 1 |
---|
145 | graphicatt "xylimits=1250,1500,0,3" |
---|
146 | plot2d dab01 (n/8192)*250+1250 val 1 "blue cpts notit nstat" |
---|
147 | plot2d dab11 (n/8192)*250+1250 val 1 "same red cpts notit nstat" |
---|
148 | settitle "Mean First DAB $source $date Ch0 blue Ch1 red" |
---|
149 | |
---|
150 | newwin 1 1 |
---|
151 | graphicatt "xylimits=1250,1500,0,3" |
---|
152 | plot2d dab02 (n/8192)*250+1250 val 1 "blue cpts notit nstat" |
---|
153 | plot2d dab12 (n/8192)*250+1250 val 1 "same red cpts notit nstat" |
---|
154 | settitle "Mean Second DAB $source $date Ch0 blue Ch1 red" |
---|
155 | |
---|
156 | ## Save results in ppf files |
---|
157 | del mtxbase mtxdab1 mtxdab2 |
---|
158 | |
---|
159 | c++exec TMatrix<r_4> mtxbase(2,base0.NElts()); mtxbase.Row(0)=base0; mtxbase.Row(1)=base1; KeepObj(mtxbase); \ |
---|
160 | TMatrix<r_4> mtxdab1(2,dab01.NElts()); mtxdab1.Row(0)=dab01; mtxdab1.Row(1)=dab11; KeepObj(mtxdab1); \ |
---|
161 | TMatrix<r_4> mtxdab2(2,dab02.NElts()); mtxdab2.Row(0)=dab02; mtxdab2.Row(1)=dab12; KeepObj(mtxdab2); |
---|
162 | |
---|
163 | saveppf mtxbase ${toppath}/${source}/${date}${srclower}/baseline_${date}_${srclower}_${mode}.ppf |
---|
164 | saveppf mtxdab1 ${toppath}/${source}/${date}${srclower}/dab1_${date}_${srclower}_${mode}.ppf |
---|
165 | saveppf mtxdab2 ${toppath}/${source}/${date}${srclower}/dab2_${date}_${srclower}_${mode}.ppf |
---|