Changeset 3769 in Sophya for trunk/Cosmo/RadioBeam/mdish.cc

Timestamp:

May 7, 2010, 6:44:43 PM (15 years ago)

Author:

ansari

Message:

Corrections/amelioration du programme de calcul de la sensibilite interfero, Reza 07/05/2010

File:

• trunk/Cosmo/RadioBeam/mdish.cc (modified) (8 diffs)

trunk/Cosmo/RadioBeam/mdish.cc

@@ -108 +108 @@
   : nx(0),ny(0),xmin(0),xmax(0),ymin(0),ymax(0),sumw0(0.)
 {
+  ixb0 = jyb0 = 0;
 }
 QHis2D::QHis2D(r_8 xMin,r_8 xMax,int_4 nxb,r_8 yMin,r_8 yMax,int_4 nyb)
@@ -123 +124 @@
   sa_size_t sz[5];  sz[0]=nx;  sz[1]=ny;
   aw.ReSize(2,sz);
+  SetZeroBin();
   sumw0=0.;
   return;
@@ -131 +133 @@
   sa_size_t jy = (sa_size_t)((y-ymin)/dyb);
   if ((ix<0)||(ix>=nx)||(jy<0)||(jy>=ny))  return 0.;
-  double rw = ((ix==0)&&(jy==0)) ? w : 0.;
+  double rw = ((ix==ixb0)&&(jy==jyb0)) ? w : 0.;
   sumw0 += rw; 
   if (fgfh) aw(ix,jy) += w;
   return rw;
 }
+void QHis2D::SetZeroBin(r_8 x, r_8 y)
+{
+  ixb0 = (sa_size_t)((x-xmin)/dxb);
+  jyb0 = (sa_size_t)((y-ymin)/dyb);
+}
 Histo2D QHis2D::Convert()
 {
+  int_4 imn,jmn,imx,jmx;
+  r_8 min = aw(0,0);
+  r_8 max = aw(0,0);
+  imn=jmn=imx=jmx=0;
   Histo2D h2(xmin,xmax,nx,ymin,ymax,ny);
   for(int_4 j=0; j<ny; j++) 
-    for(int_4 i=0; i<nx; i++)  h2(i,j) = aw(i,j);
+    for(int_4 i=0; i<nx; i++) {
+      h2(i,j) = aw(i,j);
+      if (aw(i,j)>max) {
+        imx=i;  jmx=j;  max=aw(i,j);
+      }
+      if (aw(i,j)<min) {
+        imn=i;  jmn=j;  min=aw(i,j);
+      }
+    }
+  cout << "QHis2D::Convert()/Info: Nx,Ny=" << nx << "," << ny << " SumW=" << sumw0 
+       << "\n ... Max:" << imx << "," << jmx << " ->" << max 
+       << " @" << imx*dxb+xmin << "," << jmx*dyb+ymin
+       << "\n ...Min:" << imn << "," << jmn << " ->" << min 
+       << " @" << imn*dxb+xmin << "," << jmn*dyb+ymin << endl;
   return h2;
 }
@@ -159 +183 @@
   cout << " MultiDish::GetResponse() - NDishes=" << dishes_.size() << " nx=" << nx_ << " ny=" << ny_ << endl;
   double kmx = 1.2*DeuxPI*dmax_/lambda_;
-  /*
-  h2w_= Histo2D(0.,kmx,nx_,0.,kmx,ny_);
-  h2cnt_= Histo2D(0.,kmx,nx_,0.,kmx,ny_);
-  h2w_.Zero();
-  h2cnt_.Zero();
-  */
-  h2w_.Define(0.,kmx,nx_,0.,kmx,ny_);
+  double dkmx = kmx/(double)nx_;
+  double dkmy = kmx/(double)ny_;
+  double kmxx = ((double)nx_+0.5)*dkmx;
+  double kmxy = ((double)ny_+0.5)*dkmy;
+  h2w_.Define(-kmxx,kmxx,2*nx_+1,-kmxy,kmxy,2*ny_+1);
+  h2w_.SetZeroBin(0.,0.);
 
   double dold = dishes_[0].D/lambda_;
@@ -184 +207 @@
       sumw += CumulResp(rd, (double)kt*dtet, (double)jp*dphi);
 
-  double kx0 = DeuxPI*fabs(dishes_[1].X-dishes_[0].X)/lambda_;
-  double ky0 = DeuxPI*fabs(dishes_[1].Y-dishes_[0].Y)/lambda_;
-  int_4 ib, jb;
+  double kx1 = DeuxPI*(dishes_[0].DiameterX())/lambda_;
+  double ky1 = DeuxPI*(dishes_[0].DiameterY())/lambda_;
+  int_4 ib,jb;
   //  h2w_ /= h2cnt_; 
   Histo2D h2 = h2w_.Convert();
-  h2.FindBin(kx0, ky0, ib, jb);
-  cout << " ---- MultiDish::GetResponse() VMin=" << h2.VMin() << " VMax= " << h2.VMax() 
-       << " h(0,0)=" << h2(0,0) << " h(" << ib <<"," << jb << ")=" << h2(ib,jb) <<endl;
+  h2.FindBin(kx1, ky1, ib, jb);
+  if ((kx1<0)||(ky1<0)||(kx1>=h2.NBinX())||(ky1>=h2.NBinY())) {
+    cout << " MultiDish::GetResponse[1]/ERROR kx1,ky1=" << kx1 <<","<< ky1 << " --> ib,jb=" << ib <<","<< jb << endl;
+    ib=jb=0;
+  }
+  double vmax=h2.VMax();
+  cout << " MultiDish::GetResponse[1] VMin=" << h2.VMin() << " VMax= " << vmax  
+       << " h(0,0)=" << h2(0,0) << " kx1,ky1->h(" << ib <<"," << jb << ")=" << h2(ib,jb) <<endl;
   //  double fnorm=sqrt((double)dishes_.size())/h2.VMax();
   double fnorm=1.;
-  if (h2.VMax() > sumw) {
+  if (vmax > sumw) {
     fnorm=(double)dishes_.size()/h2.VMax();
-    cout << " ---- MultiDish::GetResponse() NDishes=" << dishes_.size() << " sumw=" << sumw 
-         << " Renormalizing x NDishes/sumw  " << fnorm << endl;
+    cout << " MultiDish::GetResponse[2]/Warning h2.VMax()=" << vmax << " >  sumw=" << sumw << endl;  
+    cout << "   ... NDishes=" << dishes_.size() << " sumw=" << sumw 
+         << " Renormalizing x NDishes/VMax  " << fnorm << endl;
   }
   else {
-    fnorm=(double)dishes_.size()/h2.VMax();
-    cout << " ---- MultiDish::GetResponse() NDishes=" << dishes_.size() << " VMax=" << h2.VMax()  
-         << " Renormalizing x NDishes/h2.VMax()   " << fnorm << endl;
+    fnorm=(double)dishes_.size()/sumw;
+    cout << " MultiDish::GetResponse[3] NDishes=" << dishes_.size() << " sumw=" << sumw  
+         << " Renormalizing x NDishes/sumw   " << fnorm << endl;
   }
   h2 *= fnorm;
-  cout << " ---- MultiDish::GetResponse() APRES VMin=" << h2.VMin() << " VMax= " << h2.VMax() << " h(0,0)=" 
+  cout << " ---- MultiDish::GetResponse/[4] APRES VMin=" << h2.VMin() << " VMax= " << h2.VMax() << " h(0,0)=" 
        << h2(0,0) << endl;
   return h2;
 }
 
-/*
-double MultiDish::AddToHisto(double kx0, double ky0, double x, double y, double w, bool fgfh)
-{
-  double xxp = kx0+x;
-  double yyp = ky0+y;
-  double sumw=0.;
-  int_4 ib, jb;
-  h2w_.FindBin(xxp, yyp, ib, jb);
-  if ((ib==0)&&(jb==0))  sumw+=w;
-  if (fgfh) { 
-    h2w_.Add(xxp, yyp, w);
-    h2cnt_.Add(xxp, yyp, 1.);
-  }
-  double xxm=kx0-x;
-  double yym=ky0-y;
-  if (xxm>0.)  {
-    h2w_.FindBin(xxm, yyp, ib, jb);
-    if ((ib==0)&&(jb==0))  sumw+=w;
-    if (fgfh) {
-      h2w_.Add(xxm, yyp, w);
-      h2cnt_.Add(xxm, yyp, 1.);
-    }
-    if (yym>0.) { 
-      h2w_.FindBin(xxm, yym, ib, jb);
-      if ((ib==0)&&(jb==0))  sumw+=w;
-      if (fgfh) { 
-        h2w_.Add(xxm, yym, w);
-        h2cnt_.Add(xxm, yym, 1.);
-      }
-    }
-  }
-  if (yym>0.)  { 
-    h2w_.FindBin(xxp, yym, ib, jb);
-    if ((ib==0)&&(jb==0))  sumw+=w;
-    if (fgfh) {
-      h2w_.Add(xxp, yym, w);
-      h2cnt_.Add(xxp, yym, 1.);
-    }
-  }
-  return sumw; 
-}
-*/
+Histo2D MultiDish::PosDist(int nx, int ny, double dmax)
+{
+  if (dmax<1e-3)  dmax=nx*dishes_[0].Diameter();
+  double dd = dmax/nx/2.;
+  Histo2D hpos(-dd,dmax+dd,nx+1,-dd,dmax+dd,ny+1);
+  for(size_t i=0; i<NbDishes(); i++) {
+    hpos.Add(dishes_[i].X, dishes_[i].Y);
+  }
+  return hpos;
+}
 
 double MultiDish::AddToHisto(double kx0, double ky0, double x, double y, double w, bool fgfh)
@@ -261 +258 @@
   double xxm=kx0-x;
   double yym=ky0-y;
-  if (xxm>0.)  {
-    sumw += h2w_.Add(xxm, yyp, w, fgfh);
-    if (yym>0.)  sumw += h2w_.Add(xxm, yym, w, fgfh);
-  }
-  if (yym>0.)  sumw += h2w_.Add(xxp, yym, w, fgfh);
+  //  if (xxm>0.)  {
+  sumw += h2w_.Add(xxm, yyp, w, fgfh);
+  // if (yym>0.)  
+  sumw += h2w_.Add(xxm, yym, w, fgfh);
+  //  }
+  // if (yym>0.)  
+  sumw += h2w_.Add(xxp, yym, w, fgfh);
   return sumw; 
 }
@@ -275 +274 @@
   double dx = h2w_.WBinX()/5;
   double dy = h2w_.WBinY()/5;
-  int nbx = DeuxPI*rd.Dx()/dx;
-  int nby = DeuxPI*rd.Dy()/dy;
+  int nbx = DeuxPI*rd.Dx()/dx+1;
+  int nby = DeuxPI*rd.Dy()/dy+1;
   dx = DeuxPI*rd.Dx()/(double)nbx;
   dy = DeuxPI*rd.Dy()/(double)nby;
@@ -287 +286 @@
 
   for(size_t i=0; i<dishes_.size(); i++) {
-    for(size_t j=i; j<dishes_.size(); j++) {
-      double kx0 = DeuxPI*fabs(dishes_[i].X-dishes_[j].X)/lambda_;
-      double ky0 = DeuxPI*fabs(dishes_[i].Y-dishes_[j].Y)/lambda_;
+    for(size_t j=0; j<dishes_.size(); j++) {
+      double kx0 = DeuxPI*(dishes_[i].X-dishes_[j].X)/lambda_;
+      double ky0 = DeuxPI*(dishes_[i].Y-dishes_[j].Y)/lambda_;
       rot.Do(kx0, ky0);
-      if (kx0<0) kx0=-kx0;
-      if (ky0<0) ky0=-ky0;
+      //    if (kx0<0) kx0=-kx0;
+      //    if (ky0<0) ky0=-ky0;
       bool fgfh= (!fgnoauto_ || (dishes_[i].ReflectorId()!=dishes_[j].ReflectorId()));
       for(int ix=0; ix<nbx; ix++) 
         for(int jy=0; jy<nby; jy++) { 
           double x = ix*dx;  
-          double y = jy*dy;  
-          sumw += AddToHisto(kx0, ky0, x, y, rd(x,y), fgfh);
+          double y = jy*dy;
+          if ((ix>0)&&(jy>0)) {
+            sumw += AddToHisto(kx0, ky0, x, y, rd(x,y), fgfh);
+          }
+          else {
+            if ((ix==0)&&(jy==0)) 
+              sumw += h2w_.Add(kx0, ky0, rd(0.,0.), fgfh);
+            else {
+              double w = rd(x,y);
+              if (ix==0) {
+                sumw += h2w_.Add(kx0, ky0+y, w, fgfh);
+                sumw += h2w_.Add(kx0, ky0-y, w, fgfh);
+              }
+              else {
+                sumw += h2w_.Add(kx0+x, ky0, w, fgfh);
+                sumw += h2w_.Add(kx0-x, ky0, w, fgfh);
+              }
+            }
+            //   
+          }
         }
-    }
     //    if (i%10==0) 
     //      cout << " MultiDish::CumulResp() done i=" << i << " / imax=" << dishes_.size() 
     //     << " theta=" << theta << " phi=" << phi << endl;
+    }
   }
   return sumw;