source: Sophya/trunk/Cosmo/RadioBeam/mdish.cc@ 3973

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

Corrections diverses: choix lobe gaussien/triangle et specif DishDiameter au lieu de DoL ds applobe/calcpk2, possibilite application lobe freq.independante ds applobe, Reza 18/04/2011
File size: 13.7 KB
Line 
1// Classes to compute 2D
2// R. Ansari - Nov 2008, May 2010
3
4#include "mdish.h"
5
6
7//--------------------------------------------------
8// -- Four2DResponse class
9//--------------------------------------------------
10// Constructor
11Four2DResponse::Four2DResponse(int typ, double dx, double dy, double lambda)
12 : typ_(typ), dx_((dx>1.e-3)?dx:1.), dy_((dy>1.e-3)?dy:1.)
13{
14 setLambdaRef(lambda);
15 setLambda(lambda);
16}
17
18// Return the response for the wave vecteor (kx,ky)
19double Four2DResponse::Value(double kx, double ky)
20{
21 kx *= lambda_ratio_;
22 ky *= lambda_ratio_;
23 double wk,wkx,wky;
24 switch (typ_)
25 {
26 case 1: // Reponse gaussienne parabole diametre D exp[ -2 (lambda k_g / D )^2 ]
27 wk = sqrt(kx*kx+ky*ky)/dx_;
28 wk = 2*wk*wk;
29 return exp(-wk);
30 break;
31 case 2: // Reponse parabole diametre D Triangle <= kmax= 2 pi D / lambda
32 wk = sqrt(kx*kx+ky*ky)/dx_/2./M_PI;
33 return ( (wk<1.)?(1.-wk):0.);
34 break;
35 case 22: // Reponse parabole diametre D Triangle <= kmax= 2 pi D / lambda + trou au centre
36 wk = sqrt(kx*kx+ky*ky)/dx_/2./M_PI;
37 if (wk<0.025) return 39.*wk;
38 else if (wk<1.) return (1.-wk);
39 else return 0.;
40 break;
41 case 3: // Reponse rectangle Dx x Dy Triangle (|kx|,|k_y|) <= (2 pi Dx / lambda, 2 pi Dx / lambda)
42 wkx = fabs(kx)/2./M_PI/dx_;
43 wky = fabs(ky)/2./M_PI/dy_;
44 return ( ((wkx<1.)&&(wky<1.))?((1.-wkx)*(1-wky)):0.);
45 break;
46 default:
47 return 1.;
48 }
49}
50// Return a vector representing the power spectrum (for checking)
51Histo2D Four2DResponse::GetResponse(int nx, int ny)
52{
53 double kxmx = 1.2*DeuxPI*dx_;
54 double kymx = 1.2*DeuxPI*dy_;
55 if (typ_<3) kymx=kxmx;
56 Histo2D h2(-kxmx,kxmx,nx,-kymx,kymx,ny);
57
58 double xbc,ybc;
59 for(int_4 j=0; j<h2.NBinY(); j++)
60 for(int_4 i=0; i<h2.NBinX(); i++) {
61 h2.BinCenter(i,j,xbc,ybc);
62 h2(i,j) = Value(xbc,ybc);
63 }
64 return h2;
65}
66
67HProf Four2DResponse::GetProjNoiseLevel(int nbin, bool fgnorm1)
68{
69 Histo2D h2w = GetResponse(2*nbin, 2*nbin);
70 r_8 vmin=h2w.VMin();
71 r_8 vmax=h2w.VMax();
72 double seuil=vmax/10000.;
73 if (seuil<1.e-6) seuil=1.e-6;
74 r_8 facnorm=((fgnorm1)?vmax:1.);
75 cout << " Four2DResponse::GetProjNoiseLevel Min,Max=" << vmin << " , " << vmax
76 << " facnorm=" << facnorm << " seuil=" << seuil << endl;
77 double kmax=2.*M_PI*D();
78 HProf hp(0,kmax,nbin);
79 double x,y;
80 for(sa_size_t j=0; j<h2w.NBinY(); j++) {
81 for(sa_size_t i=0; i<h2w.NBinX(); i++) {
82 h2w.BinCenter(i,j,x,y);
83 double yw=h2w(i,j);
84 if (yw<seuil) continue;
85 hp.Add(sqrt(x*x+y*y),facnorm/yw);
86 }
87 }
88 return hp;
89}
90
91HProf Four2DResponse::GetProjResponse(int nbin, bool fgnorm1)
92{
93 Histo2D h2w = GetResponse(2*nbin, 2*nbin);
94 r_8 vmin=h2w.VMin();
95 r_8 vmax=h2w.VMax();
96 r_8 facnorm=((fgnorm1)?(1./vmax):1.);
97 cout << " Four2DResponse::GetProjResponse Min,Max=" << vmin << " , " << vmax
98 << " facnorm=" << facnorm << endl;
99 double kmax=2.*M_PI*D();
100 HProf hp(0,kmax,nbin);
101 double x,y;
102 for(sa_size_t j=0; j<h2w.NBinY(); j++) {
103 for(sa_size_t i=0; i<h2w.NBinX(); i++) {
104 h2w.BinCenter(i,j,x,y);
105 hp.Add(sqrt(x*x+y*y),h2w(i,j)*facnorm);
106 }
107 }
108 return hp;
109}
110
111//---------------------------------------------------------------
112// -- Four2DRespTable : Reponse tabulee instrumentale ds le plan k_x,k_y (angles theta,phi)
113//---------------------------------------------------------------
114Four2DRespTable::Four2DRespTable()
115 : Four2DResponse(0,1.,1.)
116{
117}
118
119Four2DRespTable::Four2DRespTable(Histo2D const & hrep, double d, double lambda)
120 : Four2DResponse(0,d,d,lambda) , hrep_(hrep)
121{
122}
123
124double Four2DRespTable::Value(double kx, double ky)
125{
126 kx *= lambda_ratio_;
127 ky *= lambda_ratio_;
128 int_4 i,j;
129 if ( (kx<=hrep_.XMin())||(kx>=hrep_.XMax()) ||
130 (ky<=hrep_.YMin())||(ky>=hrep_.YMax()) ) return 0.;
131 hrep_.FindBin(kx, ky, i, j);
132 return hrep_(i, j);
133}
134
135double Four2DRespTable::renormalize(double max)
136{
137 double cmx = hrep_.VMax();
138 hrep_ *= (max/cmx);
139 return cmx;
140}
141
142void Four2DRespTable::writeToPPF(string flnm)
143{
144 DVList dvinfo;
145 dvinfo["DoL"] = dx_;
146 dvinfo["LambdaRef"] = lambdaref_;
147 dvinfo["Lambda"] = lambda_;
148 POutPersist po(flnm);
149 po << hrep_;
150 po << dvinfo;
151}
152
153void Four2DRespTable::readFromPPF(string flnm)
154{
155 PInPersist pin(flnm);
156 DVList dvinfo;
157 pin >> hrep_;
158 pin >> dvinfo;
159 dx_ = dy_ = dvinfo["DoL"];
160 setLambdaRef((double)dvinfo["LambdaRef"]);
161 setLambda((double)dvinfo["Lambda"]);
162}
163
164
165
166//---------------------------------------------------------------
167// -- Four2DRespRatio : rapport de la reponse entre deux objets Four2DResponse
168//---------------------------------------------------------------
169Four2DRespRatio::Four2DRespRatio(Four2DResponse& a, Four2DResponse& b, double divzthr)
170 : Four2DResponse(0, a.D(), a.D()), a_(a), b_(b), divzthr_(divzthr)
171{
172}
173
174double Four2DRespRatio::Value(double kx, double ky)
175{
176 double ra = a_.Value(kx,ky);
177 double rb = b_.Value(kx,ky);
178 if (ra<rb) {
179 if (rb>1.e-39) return(ra/rb);
180 else return 0.;
181 }
182 if (rb<divzthr_) rb=divzthr_;
183 return (ra/rb);
184}
185
186//---------------------------------------------------------------
187//--- Classe simple pour le calcul de rotation
188class Rotation {
189public:
190 Rotation(double tet, double phi)
191 {
192// (Teta,Phi) = Direction de visee
193// Les angles d'Euler correspondants sont Teta, Phi+Pi/2
194// Le Pi/2 vient que les rotations d'euler se font dans l'ordre
195// Autour de oZ d'angle Phi, autour de oN (nouvel axe X) d'angle Teta
196// Autour du nouvel axe Z (x3) d'angle Psi (Psi=0 -> cp=1, sp=0.;
197 double ct = cos(tet);
198 double st = sin(tet);
199 // Le Pi/2 echange les axes X<>Y pour theta=phi=0 !
200 // double cf = cos(phi+M_PI/2);
201 // double sf = sin(phi+M_PI/2);
202 double cf = cos(phi);
203 double sf = sin(phi);
204 double cp = 1.; // cos((double)pO);
205 double sp = 0.; // sin((double)pO);
206 RE[0][0] = cf*cp-sf*ct*sp; RE[0][1] = sf*cp+cf*ct*sp; RE[0][2] = st*sp;
207 RE[1][0] = -cf*sp-sf*ct*cp; RE[1][1] = -sf*sp+cf*ct*cp; RE[1][2] = st*cp;
208 RE[2][0] = sf*st; RE[2][1] = -cf*st; RE[2][2] = ct;
209 }
210 inline void Do(double& x, double& y)
211 {
212 double xx=x;
213 double yy=y;
214 x = RE[0][0]*xx+RE[0][1]*yy;
215 y = RE[1][0]*xx+RE[1][1]*yy;
216 }
217 double RE[3][3];
218};
219
220
221//----------------------------------------------------------------------
222// -- Pour calculer la reponse ds le plan kx,ky d'un system MultiDish
223//----------------------------------------------------------------------
224MultiDish::MultiDish(double lambda, double dmax, vector<Dish>& dishes, bool fgnoauto)
225 : lambda_(lambda), dmax_(dmax), dishes_(dishes), fgnoauto_(fgnoauto)
226{
227 SetThetaPhiRange();
228 SetRespHisNBins();
229 SetBeamNSamples();
230 SetPrtLevel();
231 fgcomputedone_=false;
232 mcnt_=0;
233}
234
235void MultiDish::ComputeResponse()
236{
237 cout << " MultiDish::ComputeResponse() - NDishes=" << dishes_.size() << " nx=" << nx_ << " ny=" << ny_ << endl;
238 double kmx = 1.2*DeuxPI*dmax_/lambda_;
239 double dkmx = kmx/(double)nx_;
240 double dkmy = kmx/(double)ny_;
241 double kmxx = ((double)nx_+0.5)*dkmx;
242 double kmxy = ((double)ny_+0.5)*dkmy;
243 h2w_.Define(-kmxx,kmxx,2*nx_+1,-kmxy,kmxy,2*ny_+1);
244 h2w_.SetZeroBin(0.,0.);
245 kmax_=kmx;
246
247 double dold = dishes_[0].Diameter()/lambda_;
248 double dolx = dishes_[0].DiameterX()/lambda_;
249 double doly = dishes_[0].DiameterY()/lambda_;
250
251 Four2DResponse rd(2, dold, dold);
252 Four2DResponse rdr(3, dolx, doly);
253
254 if (!dishes_[0].isCircular()) rd = rdr;
255
256 double dtet = thetamax_/(double)ntet_;
257 double dphi = phimax_/(double)nphi_;
258 cout << " MultiDish::ComputeResponse() - ThetaMax=" << thetamax_ << " NTheta=" << ntet_
259 << " PhiMax=" << phimax_ << " NPhi=" << nphi_ << endl;
260
261 double sumw = 0.;
262 for(int kt=0; kt<ntet_; kt++) {
263 double theta=(double)kt*dtet;
264 for(int jp=0; jp<nphi_; jp++) {
265 double phi=(double)jp*dphi;
266 sumw += CumulResp(rd, theta, phi);
267 if (theta<1.e-9) continue;
268 sumw += CumulResp(rd, theta, -phi);
269 sumw += CumulResp(rd, theta, phi+M_PI);
270 sumw += CumulResp(rd, theta, -(phi+M_PI));
271 }
272 if (prtlev_>0)
273 cout << " MultiDish::ComputeResponse() done ktheta=" << kt << " / MaxNTheta= "
274 << ntet_ << endl;
275 }
276 r_8 rmin,rmax;
277 h2w_.GetMinMax(rmin,rmax);
278 cout << " MultiDish::ComputeResponse() finished : Rep_min,max=" << rmin << "," << rmax << " sumW0="
279 << sumw << " ?=" << h2w_.SumWBinZero() << endl;
280 fgcomputedone_=true;
281 return;
282}
283
284Histo2D MultiDish::GetResponse()
285{
286 if (!fgcomputedone_) ComputeResponse();
287
288 double kx1 = DeuxPI*(dishes_[0].DiameterX())/lambda_;
289 double ky1 = DeuxPI*(dishes_[0].DiameterY())/lambda_;
290 int_4 ib,jb;
291 // h2w_ /= h2cnt_;
292 Histo2D h2 = h2w_.Convert();
293 h2.FindBin(kx1, ky1, ib, jb);
294 if ((kx1<0)||(ky1<0)||(kx1>=h2.NBinX())||(ky1>=h2.NBinY())) {
295 cout << " MultiDish::GetResponse[1]/ERROR kx1,ky1=" << kx1 <<","<< ky1 << " --> ib,jb=" << ib <<","<< jb << endl;
296 ib=jb=0;
297 }
298 double sumw=h2w_.SumWBinZero();
299 double vmax=h2.VMax();
300 cout << " MultiDish::GetResponse[1] VMin=" << h2.VMin() << " VMax= " << vmax
301 << " h(0,0)=" << h2(0,0) << " kx1,ky1->h(" << ib <<"," << jb << ")=" << h2(ib,jb) <<endl;
302 // double fnorm=sqrt((double)dishes_.size())/h2.VMax();
303 double fnorm=1.;
304 if (vmax > sumw) {
305 fnorm=(double)dishes_.size()/h2.VMax();
306 cout << " MultiDish::GetResponse[2]/Warning h2.VMax()=" << vmax << " > sumw=" << sumw << endl;
307 cout << " ... NDishes=" << dishes_.size() << " sumw=" << sumw
308 << " Renormalizing x NDishes/VMax " << fnorm << endl;
309 }
310 else {
311 fnorm=(double)dishes_.size()/sumw;
312 cout << " MultiDish::GetResponse[3] NDishes=" << dishes_.size() << " sumw=" << sumw
313 << " Renormalizing x NDishes/sumw " << fnorm << endl;
314 }
315 h2 *= fnorm;
316 cout << " ---- MultiDish::GetResponse/[4] APRES VMin=" << h2.VMin() << " VMax= " << h2.VMax() << " h(0,0)="
317 << h2(0,0) << endl;
318 return h2;
319}
320
321HProf MultiDish::GetProjNoiseLevel(int nbin, bool fgnorm1)
322{
323 r_8 vmin,vmax;
324 h2w_.GetMinMax(vmin,vmax);
325 double seuil=vmax/10000.;
326 if (seuil<1.e-6) seuil=1.e-6;
327 r_8 facnorm=((fgnorm1)?vmax:1.);
328 cout << " MultiDish::GetProjNoiseLevel Min,Max=" << vmin << " , " << vmax
329 << " facnorm=" << facnorm << " seuil=" << seuil << endl;
330 HProf hp(0,kmax_,nbin);
331 for(sa_size_t j=0; j<h2w_.NBinY(); j++) {
332 double y=h2w_.Y(j);
333 for(sa_size_t i=0; i<h2w_.NBinX(); i++) {
334 double x=h2w_.X(i);
335 double yw=h2w_(i,j);
336 if (yw<seuil) continue;
337 hp.Add(sqrt(x*x+y*y),facnorm/yw);
338 }
339 }
340 return hp;
341}
342
343HProf MultiDish::GetProjResponse(int nbin, bool fgnorm1)
344{
345 r_8 vmin,vmax;
346 h2w_.GetMinMax(vmin,vmax);
347 r_8 facnorm=((fgnorm1)?(1./vmax):1.);
348 cout << " MultiDish::GetProjResponse Min,Max=" << vmin << " , " << vmax
349 << " facnorm=" << facnorm << endl;
350 HProf hp(0,kmax_,nbin);
351 for(sa_size_t j=0; j<h2w_.NBinY(); j++) {
352 double y=h2w_.Y(j);
353 for(sa_size_t i=0; i<h2w_.NBinX(); i++) {
354 double x=h2w_.X(i);
355 hp.Add(sqrt(x*x+y*y),h2w_(i,j)*facnorm);
356 }
357 }
358 return hp;
359}
360
361
362Histo2D MultiDish::PosDist(int nx, int ny, double dmax)
363{
364 if (dmax<1e-3) dmax=nx*dishes_[0].Diameter();
365 double dd = dmax/nx/2.;
366 Histo2D hpos(-dd,dmax+dd,nx+1,-dd,dmax+dd,ny+1);
367 for(size_t i=0; i<NbDishes(); i++) {
368 hpos.Add(dishes_[i].X, dishes_[i].Y);
369 }
370 return hpos;
371}
372
373double MultiDish::AddToHisto(double kx0, double ky0, double x, double y, double w, bool fgfh)
374{
375 double xxp = kx0+x;
376 double yyp = ky0+y;
377 double sumw=0.;
378 sumw += h2w_.Add(xxp, yyp, w, fgfh);
379 double xxm=kx0-x;
380 double yym=ky0-y;
381 // if (xxm>0.) {
382 sumw += h2w_.Add(xxm, yyp, w, fgfh);
383 // if (yym>0.)
384 sumw += h2w_.Add(xxm, yym, w, fgfh);
385 // }
386 // if (yym>0.)
387 sumw += h2w_.Add(xxp, yym, w, fgfh);
388 return sumw;
389}
390
391double MultiDish::CumulResp(Four2DResponse& rd, double theta, double phi)
392{
393 // cout << " MultiDish::CumulResp() theta=" << theta << " phi=" << phi << endl;
394 /*
395 double dx = h2w_.WBinX()/5;
396 double dy = h2w_.WBinY()/5;
397 int nbx = DeuxPI*rd.Dx()/dx+1;
398 int nby = DeuxPI*rd.Dy()/dy+1;
399 */
400 double dx,dy;
401 int nbx=beamnx_;
402 int nby=beamny_;
403 dx = DeuxPI*rd.Dx()/(double)nbx;
404 dy = DeuxPI*rd.Dy()/(double)nby;
405 if (mcnt_==0)
406 cout << " CumulResp() nbx=" << nbx << " nby=" << nby << " dx=" << dx << " dy=" << dy << endl;
407 mcnt_++;
408
409 double sumw = 0.;
410 Rotation rot(theta, phi);
411
412 for(size_t i=0; i<dishes_.size(); i++) {
413 for(size_t j=i; j<dishes_.size(); j++) {
414 double kx0 = DeuxPI*(dishes_[i].X-dishes_[j].X)/lambda_;
415 double ky0 = DeuxPI*(dishes_[i].Y-dishes_[j].Y)/lambda_;
416 double pgain=dishes_[i].Gain()*dishes_[j].Gain();
417 rot.Do(kx0, ky0);
418 // if (kx0<0) kx0=-kx0;
419 // if (ky0<0) ky0=-ky0;
420 bool fgfh= (!fgnoauto_ || (dishes_[i].ReflectorId()!=dishes_[j].ReflectorId()));
421 for(int ix=0; ix<nbx; ix++)
422 for(int jy=0; jy<nby; jy++) {
423 double x = ix*dx;
424 double y = jy*dy;
425 if ((ix>0)&&(jy>0)) {
426 sumw += AddToHisto(kx0, ky0, x, y, rd(x,y)*pgain, fgfh);
427 if (j!=i) sumw += AddToHisto(-kx0, -ky0, x, y, rd(x,y)*pgain, fgfh);
428 }
429 else {
430 if ((ix==0)&&(jy==0)) {
431 sumw += h2w_.Add(kx0, ky0, rd(0.,0.)*pgain, fgfh);
432 if (j!=i) sumw += h2w_.Add(-kx0, -ky0, rd(0.,0.)*pgain, fgfh);
433 }
434 else {
435 double w = rd(x,y)*pgain;
436 if (ix==0) {
437 sumw += h2w_.Add(kx0, ky0+y, w, fgfh);
438 sumw += h2w_.Add(kx0, ky0-y, w, fgfh);
439 if (j!=i) {
440 sumw += h2w_.Add(-kx0, -ky0+y, w, fgfh);
441 sumw += h2w_.Add(-kx0, -ky0-y, w, fgfh);
442 }
443 }
444 else {
445 sumw += h2w_.Add(kx0+x, ky0, w, fgfh);
446 sumw += h2w_.Add(kx0-x, ky0, w, fgfh);
447 if (j!=i) {
448 sumw += h2w_.Add(-kx0+x, -ky0, w, fgfh);
449 sumw += h2w_.Add(-kx0-x, -ky0, w, fgfh);
450 }
451 }
452 }
453 //
454 }
455 }
456 // if (i%10==0)
457 // cout << " MultiDish::CumulResp() done i=" << i << " / imax=" << dishes_.size()
458 // << " theta=" << theta << " phi=" << phi << endl;
459 }
460 }
461 return sumw;
462}
463
Note: See TracBrowser for help on using the repository browser.