source: Sophya/trunk/Cosmo/RadioBeam/fgndsub.cc@ 4060

Last change on this file since 4060 was 4022, checked in by ansari, 14 years ago

modifs pour pouvoir imposer la moyenne en temp des plans X,Y des cubes lors de l'extraction du signal HI, Reza 28/9/2011
File size: 11.1 KB
Line 
1/* ------------------------ Projet BAORadio --------------------
2 Classe ForegroundCleaner
3 R. Ansari , C. Magneville - Juin 2010
4--------------------------------------------------------------- */
5
6#include "fgndsub.h"
7#include "lobe.h"
8
9#include "cubedef.h"
10#include "matharr.h"
11#include "poly.h"
12
13#include "ctimer.h"
14
15/* --Methode-- */
16PowerLawChecker::PowerLawChecker(TArray< TF >& skycube)
17 : skycube_(skycube)
18{
19}
20/* --Methode-- */
21void PowerLawChecker::CheckXYMean()
22{
23 // double freq0 : Frequence premier index en k (MHz)
24 // double dfreq : // largeur en frequence de chaque plan (Mhz)
25 double freq0_ = Freq0MHz;
26 double dfreq_ = FreqSizeMHz/(double)NFreq;
27
28 double tempfirst,templast;
29 r_8 s1, sx, sx2, sy, sxy;
30 s1=sx=sx2=sy=sxy=0.;
31 double lnf0=0.;
32 for(sa_size_t k=0; k<skycube_.SizeZ(); k++) {
33 double lnf=log((double)k*dfreq_+freq0_);
34 double ttot = Mean(skycube_(Range::all(), Range::all(), Range(k)));
35 if (k==0) { tempfirst=ttot; lnf0=lnf; }
36 if (k==skycube_.SizeZ()-1) templast=ttot;
37 if (ttot < 1.e-5) continue;
38 double lntt=log(ttot);
39 s1+=1.; sx+=lnf; sx2+=(lnf*lnf);
40 sy+=lntt; sxy+=(lnf*lntt);
41 }
42 double beta = (sx*sxy-sx2*sy)/(sx*sx-s1*sx2);
43 double alpha = (s1*sxy-sx*sy)/(s1*sx2-sx*sx);
44 double T0 = exp(beta+alpha*lnf0);
45
46 cout << " PowerLawChecker::CheckMean() meanTemp(0 ...last) " << tempfirst << " ... "
47 << templast << endl;
48 bool fgnan = false;
49 if (!isfinite(alpha)||(!isfinite(beta))) {
50 cout << "ePowerLawChecker::CheckMean() Not finite alpha, beta " << endl;
51 fgnan = true;
52 }
53 cout << "PowerLawChecker::CheckMean() - T0=" << T0 << " alpha=" << alpha << "(beta=" << beta << ")" << endl;
54 return;
55}
56
57/* --Methode-- */
58ForegroundCleaner::ForegroundCleaner(Four2DResponse& arrep, Four2DResponse& tbeam, TArray< TF >& skycube, double maxratio)
59 : arrep_(arrep) , tbeam_(tbeam), skycube_(skycube), maxratio_(maxratio)
60{
61 double dxdeg = ThetaSizeDegre/(double)NTheta;
62 double dydeg = PhiSizeDegre/(double)NPhi;
63 dx_ = DegreeToRadian(dxdeg);
64 dy_ = DegreeToRadian(dydeg);
65 freq0_ = Freq0MHz;
66 dfreq_ = FreqSizeMHz/(double)NFreq;
67 cout << " ForegroundCleaner: " << " dx=" << dxdeg << " dy=" << dydeg << " degres ( dx_rad=" << dx_ << " dy_rad=" << dy_ << ")"
68 << " Freq0=" << freq0_ << " deltaFreq=" << dfreq_ << " MHz" << endl;
69 skycube.Show();
70}
71
72/* --Methode-- */
73void ForegroundCleaner::BeamCorrections()
74{
75 BeamEffect beam(arrep_);
76 beam.Correct2RefLobe(tbeam_, skycube_, dx_, dy_, freq0_, dfreq_, maxratio_);
77 cout << " ForegroundCleaner::BeamCorrections() done " << endl;
78}
79
80/* --Methode-- */
81int ForegroundCleaner::FixMeanXYTemp(double T0, double alpha)
82{
83 cout << "ForegroundCleaner::FixMeanXYTemp(T0=" << T0 << ",alpha=" << alpha << ")" << endl;
84 double lnf0=log(freq0_);
85 sa_size_t modprt=skycube_.SizeZ()/12;
86 for(sa_size_t k=0; k<skycube_.SizeZ(); k++) {
87 double lnf=log((double)k*dfreq_+freq0_);
88 double fittedtemp = T0*exp(alpha*(lnf-lnf0));
89 TArray<TF> slice = skycube_(Range::all(), Range::all(), Range(k));
90 TF deltatemp = (TF)(fittedtemp-(double)Mean(slice));
91 slice += deltatemp;
92 if (k%modprt == 0)
93 cout << "FixMeanXYTemp[k=" << k << " MeanXYTemp=" << fittedtemp-deltatemp << " -> " << fittedtemp
94 << " (DeltaTemp=" << deltatemp << ")" << endl;
95 }
96 cout << "ForegroundCleaner::FixMeanXYTemp done" << endl;
97 return 0;
98}
99
100/* --Methode-- */
101int ForegroundCleaner::CleanNegatives(TF seuil)
102{
103 sa_size_t nneg = 0.;
104 for(sa_size_t kz=0; kz<skycube_.SizeZ(); kz++)
105 for(sa_size_t ky=0; ky<skycube_.SizeY(); ky++)
106 for(sa_size_t kx=0; kx<skycube_.SizeX(); kx++)
107 if (skycube_(kx, ky, kz) < seuil) {
108 nneg++; skycube_(kx, ky, kz)=seuil;
109 }
110 cout << " ForegroundCleaner::CleanNegatives " << nneg << " sky-pixels <" << seuil << " changed to" << seuil << endl;
111 return (int)nneg;
112}
113
114/* --Methode-- */
115int ForegroundCleaner::CleanPointSources(double nsigmas)
116{
117 Timer tm("ForegroundCleaner::CleanPointSources");
118 TArray< TF > sky2d(skycube_.SizeX(), skycube_.SizeY());
119 for(sa_size_t ky=0; ky<sky2d.SizeY(); ky++)
120 for(sa_size_t kx=0; kx<sky2d.SizeX(); kx++)
121 sky2d(kx, ky) = skycube_(Range(kx), Range(ky), Range::all()).Sum();
122
123
124 double mean, sigma;
125
126
127 TArray< TF > amz(1,1,skycube_.SizeZ()), asz(1,1,skycube_.SizeZ());
128 for(sa_size_t kz=0; kz<skycube_.SizeZ(); kz++) {
129 TArray< TF > slice = skycube_(Range::all() , Range::all(), Range(kz));
130 MeanSigma(slice, mean, sigma);
131 amz(0,0,kz)=mean; asz(0,0,kz)=sigma;
132 }
133
134 MeanSigma(sky2d, mean, sigma);
135 cout << " ForegroundCleaner::CleanPointSources 2D Sky projection, mean=" << mean << " sigma=" << sigma << endl;
136
137 TF seuil = (TF)(mean+nsigmas*sigma);
138
139 sa_size_t srccnt=0;
140 for(sa_size_t ky=0; ky<skycube_.SizeY(); ky++)
141 for(sa_size_t kx=0; kx<skycube_.SizeX(); kx++) {
142 if (sky2d(kx,ky)>seuil) {
143 srccnt++;
144 skycube_(Range(kx), Range(ky), Range::all()) = amz;
145 }
146 }
147
148 cout << " Cleaned NSrc= " << srccnt << " 2D source/pixels (TotNPix=" << sky2d.Size()
149 << ")-> " << 100.*srccnt/sky2d.Size() << "% with S>" << seuil << " NSigmas=" << nsigmas << endl;
150 return (int)srccnt;
151}
152
153
154/* --Methode-- */
155TArray< TF > ForegroundCleaner::extractLSSCubeP1(TArray< TF >& synctemp, TArray< TF >& specidx)
156{
157 Timer tm("ForegroundCleaner::extractLSSCubeP1");
158// Inputs : maplss, mapsyc, freq0, dfreq
159// Outputs : synctemp, specidx (reconstructed foreground temperature and spectral index
160// Return_Array : foreground subtracted LSS signal
161 sa_size_t sz[5]; sz[0]=skycube_.SizeX(); sz[1]=skycube_.SizeY();
162 synctemp.SetSize(2, sz);
163 specidx.SetSize(2, sz);
164 TArray<r_4> omap;
165 omap.SetSize(skycube_, true);
166 Vector vlnf(skycube_.SizeZ());
167 int nprt = 0;
168 // double freq0 : Frequence premier index en k (MHz)
169 // double dfreq : // largeur en frequence de chaque plan (Mhz)
170 for(sa_size_t k=0; k<skycube_.SizeZ(); k++)
171 vlnf(k)=log((double)k*dfreq_+freq0_);
172
173 sa_size_t nbinfini=0;
174 sa_size_t nbbad=0;
175 sa_size_t imodprt=skycube_.SizeX()/6;
176 sa_size_t jmodprt=skycube_.SizeY()/6;
177 for(sa_size_t i=0; i<skycube_.SizeX(); i++)
178 for(sa_size_t j=0; j<skycube_.SizeY(); j++) {
179 r_8 s1, sx, sx2, sy, sxy;
180 s1=sx=sx2=sy=sxy=0.;
181 for(sa_size_t k=0; k<skycube_.SizeZ(); k++) {
182 double lnf = vlnf(k);
183 double ttot=(r_8)(skycube_(i,j,k));
184 if (ttot < 1.e-5) continue;
185 double lntt=log(ttot);
186 s1+=1.; sx+=lnf; sx2+=(lnf*lnf);
187 sy+=lntt; sxy+=(lnf*lntt);
188 }
189 double beta = (sx*sxy-sx2*sy)/(sx*sx-s1*sx2);
190 double alpha = (s1*sxy-sx*sy)/(s1*sx2-sx*sx);
191 double T0 = exp(beta+alpha*vlnf(0));
192
193 bool fgnan = false;
194 if (!isfinite(alpha)||(!isfinite(beta))) {
195 cout << "extractLSSCubeP1[" << i << "," << j << "]/ Not finite alpha, beta - (mapsync="
196 << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl;
197 alpha=beta=-999.;
198 fgnan = true; nbinfini++;
199 }
200 else {
201 double axp1 = beta+alpha*vlnf(0);
202 double axp2 = beta+alpha*vlnf(vlnf.Size()-1);
203
204 if ((axp1<-70.)||(axp1>70.)||(axp2<-70.)||(axp2>70.)) {
205 cout << "extractLSSCubeP1[" << i << "," << j << "] BAD alpha=" << alpha
206 << " beta=" << beta << " T0=" << T0 << " - (mapsync="
207 << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl;
208 fgnan = true; nbbad++;
209 }
210 }
211 if ((i%imodprt==0)&&(j%jmodprt==0))
212 cout << "extractLSSCubeP1[" << i << "," << j << "]: T0=" << T0 << " alpha=" << alpha
213 << " (mapsync=" << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl;
214 synctemp(i,j) = T0;
215 specidx(i,j) = alpha;
216 if (fgnan) {
217 for(sa_size_t k=0; k<skycube_.SizeZ(); k++) omap(i,j,k) = 0.;
218 }
219 else {
220 for(sa_size_t k=0; k<skycube_.SizeZ(); k++) {
221 r_4 fittedtemp = (r_4)(exp(beta+alpha*vlnf(k)));
222 omap(i,j,k) = skycube_(i,j,k)-fittedtemp;
223 }
224 }
225 }
226 cout << " ForegroundCleaner::extractLSSCubeP1() - NbNan alpha/beta=" << nbinfini << " NbBAD =" << nbbad << endl;
227 return omap;
228}
229
230/*
231static inline val_polyn2(double alpha, double beta, double gamma, double x)
232{
233 return (beta+alpha*x+gamma*x*x);
234}
235*/
236
237/* --Methode-- */
238TArray< TF > ForegroundCleaner::extractLSSCubeP2(TArray< TF >& synctemp, TArray< TF >& specidx)
239{
240 Timer tm("ForegroundCleaner::extractLSSCubeP2");
241// Inputs : maplss, mapsyc, freq0, dfreq
242// Outputs : synctemp, specidx (reconstructed foreground temperature and spectral index
243// Return_Array : foreground subtracted LSS signal
244 sa_size_t sz[5]; sz[0]=skycube_.SizeX(); sz[1]=skycube_.SizeY();
245 synctemp.SetSize(2, sz);
246 specidx.SetSize(2, sz);
247 TArray<r_4> omap;
248 omap.SetSize(skycube_, true);
249 Vector vlnf(skycube_.SizeZ());
250 Vector vlnT(skycube_.SizeZ());
251 int nprt = 0;
252 // double freq0 : Frequence premier index en k (MHz)
253 // double dfreq : // largeur en frequence de chaque plan (Mhz)
254 for(sa_size_t k=0; k<skycube_.SizeZ(); k++)
255 vlnf(k)=log((double)k*dfreq_+freq0_);
256 vlnf -= vlnf(0);
257 // cout << " DBG*extractLSSCubeP2 vlnf(0)=" << vlnf(0) << " vlnf(1)=" << vlnf(1)
258 // << "vlnf(last)=" << vlnf(vlnf.Size()-1) << endl;
259 sa_size_t nbinfini=0;
260 sa_size_t nbbad=0;
261 sa_size_t imodprt=skycube_.SizeX()/6;
262 sa_size_t jmodprt=skycube_.SizeY()/6;
263 for(sa_size_t i=0; i<skycube_.SizeX(); i++)
264 for(sa_size_t j=0; j<skycube_.SizeY(); j++) {
265 vlnT = -12.;
266 Poly polyn;
267 for(sa_size_t k=0; k<skycube_.SizeZ(); k++) {
268 double lnf = vlnf(k);
269 double ttot=(r_8)(skycube_(i,j,k));
270 if (ttot < 1.e-5) continue;
271 vlnT(k)=log(ttot);
272 }
273 polyn.Fit(vlnf,vlnT,2);
274 double beta = polyn[0];
275 double alpha = polyn[1];
276 double gamma = polyn[2];
277 double T0 = exp(polyn(vlnf(0))); // exp( val_polyn2(alpha, beta, gamma, vlnf(0)) );
278
279 bool fgnan = false;
280 if (!isfinite(alpha)||(!isfinite(beta))||(!isfinite(gamma))) {
281 cout << "extractLSSCubeP2[" << i << "," << j << "]/ Not finite alpha, beta - (mapsync="
282 << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl;
283 alpha=beta=gamma=-999.;
284 fgnan = true; nbinfini++;
285 }
286 else {
287 double axp1 = polyn(vlnf(0));
288 double axp2 = polyn(vlnf(vlnf.Size()-1));
289
290 if ((axp1<-70.)||(axp1>70.)||(axp2<-70.)||(axp2>70.)) {
291 cout << "extractLSSCubeP2[" << i << "," << j << "] BAD alpha=" << alpha
292 // << " beta=" << beta << " gamma=" << gamma << " T0=" << T0 << " axp1=" << axp1 << " axp2=" << axp2
293 << " beta=" << beta << " gamma=" << gamma << " T0=" << T0 << " - (mapsync=" << skycube_(i,j,0)
294 << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl;
295 fgnan = true; nbbad++;
296 }
297 }
298 if ((i%imodprt==0)&&(j%jmodprt==0))
299 cout << "extractLSSCubeP2[" << i << "," << j << "]: T0=" << T0 << " alpha=" << alpha << " gamma=" << gamma
300 << " (mapsync=" << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl;
301 synctemp(i,j) = T0;
302 specidx(i,j) = alpha;
303 if (fgnan) {
304 for(sa_size_t k=0; k<skycube_.SizeZ(); k++) omap(i,j,k) = 0.;
305 }
306 else {
307 for(sa_size_t k=0; k<skycube_.SizeZ(); k++) {
308 r_4 fittedtemp = (r_4)( exp(polyn(vlnf(k))) );
309 omap(i,j,k) = skycube_(i,j,k)-fittedtemp;
310 }
311 }
312 }
313 cout << " ForegroundCleaner::extractLSSCubeP2() - NbNan alpha/beta=" << nbinfini << " NbBAD =" << nbbad << endl;
314 return omap;
315}
316
Note: See TracBrowser for help on using the repository browser.