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

Timestamp:

Feb 14, 2011, 12:58:29 AM (15 years ago)

Author:

ansari

Message:

Amelioration et modifs diverses lors de la redaction du papier, Reza 13/02/2011

File:

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

trunk/Cosmo/RadioBeam/mdish.cc

@@ -65 +65 @@
 }
 
+HProf Four2DResponse::GetProjNoiseLevel(int nbin, bool fgnorm1)
+{
+  Histo2D h2w = GetResponse(2*nbin, 2*nbin);
+  r_8 vmin=h2w.VMin();
+  r_8 vmax=h2w.VMax();
+  double seuil=vmax/10000.;
+  if (seuil<1.e-6) seuil=1.e-6;
+  r_8 facnorm=((fgnorm1)?vmax:1.);
+  cout << " Four2DResponse::GetProjNoiseLevel Min,Max=" << vmin << " , " << vmax 
+       << " facnorm=" << facnorm << " seuil=" << seuil << endl;  
+  double kmax=2.*M_PI*D();
+  HProf hp(0,kmax,nbin);
+  double x,y;
+  for(sa_size_t j=0; j<h2w.NBinY(); j++) {
+    for(sa_size_t i=0; i<h2w.NBinX(); i++) {
+      h2w.BinCenter(i,j,x,y);
+      double yw=h2w(i,j);
+      if (yw<seuil) continue;
+      hp.Add(sqrt(x*x+y*y),facnorm/yw);
+    }
+  }
+  return hp;
+}
+
+HProf Four2DResponse::GetProjResponse(int nbin, bool fgnorm1)
+{
+  Histo2D h2w = GetResponse(2*nbin, 2*nbin);
+  r_8 vmin=h2w.VMin();
+  r_8 vmax=h2w.VMax();
+  r_8 facnorm=((fgnorm1)?(1./vmax):1.);
+  cout << " Four2DResponse::GetProjResponse Min,Max=" << vmin << " , " << vmax 
+       << " facnorm=" << facnorm << endl; 
+  double kmax=2.*M_PI*D();
+  HProf hp(0,kmax,nbin);
+  double x,y;
+  for(sa_size_t j=0; j<h2w.NBinY(); j++) {
+    for(sa_size_t i=0; i<h2w.NBinX(); i++) {
+      h2w.BinCenter(i,j,x,y);
+      hp.Add(sqrt(x*x+y*y),h2w(i,j)*facnorm);
+    }
+  }
+  return hp;
+}
+
 //---------------------------------------------------------------
 // -- Four2DRespTable : Reponse tabulee instrumentale ds le plan k_x,k_y (angles theta,phi) 
@@ -185 +229 @@
   SetBeamNSamples();
   SetPrtLevel();
+  fgcomputedone_=false;
   mcnt_=0;
 }
 
-Histo2D MultiDish::GetResponse()
-{
-  cout << " MultiDish::GetResponse() - NDishes=" << dishes_.size() << " nx=" << nx_ << " ny=" << ny_ << endl;
+void MultiDish::ComputeResponse()
+{
+  cout << " MultiDish::ComputeResponse() - NDishes=" << dishes_.size() << " nx=" << nx_ << " ny=" << ny_ << endl;
   double kmx = 1.2*DeuxPI*dmax_/lambda_;
   double dkmx = kmx/(double)nx_;
@@ -198 +243 @@
   h2w_.Define(-kmxx,kmxx,2*nx_+1,-kmxy,kmxy,2*ny_+1);
   h2w_.SetZeroBin(0.,0.);
+  kmax_=kmx;
 
   double dold = dishes_[0].Diameter()/lambda_;
@@ -210 +256 @@
   double dtet = thetamax_/(double)ntet_;
   double dphi = phimax_/(double)nphi_;
-  cout << " MultiDish::GetResponse() - ThetaMax=" << thetamax_ << " NTheta=" << ntet_ 
+  cout << " MultiDish::ComputeResponse() - ThetaMax=" << thetamax_ << " NTheta=" << ntet_ 
        << " PhiMax=" <<  phimax_ << " NPhi=" << nphi_ << endl;
 
@@ -224 +270 @@
       sumw += CumulResp(rd, theta, -(phi+M_PI));
     }
-    if (prtlev_>0) cout << "  MultiDish::GetResponse() done ktheta=" << kt << " / MaxNTheta= " << ntet_ << endl; 
-  }
+    if (prtlev_>0) 
+      cout << "  MultiDish::ComputeResponse() done ktheta=" << kt << " / MaxNTheta= " 
+           << ntet_ << endl; 
+  }
+  r_8 rmin,rmax;
+  h2w_.GetMinMax(rmin,rmax);
+  cout << "  MultiDish::ComputeResponse() finished : Rep_min,max=" << rmin << "," << rmax << " sumW0=" 
+       << sumw << " ?=" << h2w_.SumWBinZero() << endl;
+  fgcomputedone_=true;
+  return;
+}
+
+Histo2D MultiDish::GetResponse()
+{
+  if (!fgcomputedone_) ComputeResponse();
+
   double kx1 = DeuxPI*(dishes_[0].DiameterX())/lambda_;
   double ky1 = DeuxPI*(dishes_[0].DiameterY())/lambda_;
@@ -236 +296 @@
     ib=jb=0;
   }
+  double sumw=h2w_.SumWBinZero();
   double vmax=h2.VMax();
   cout << " MultiDish::GetResponse[1] VMin=" << h2.VMin() << " VMax= " << vmax  
@@ -257 +318 @@
   return h2;
 }
+
+HProf MultiDish::GetProjNoiseLevel(int nbin, bool fgnorm1)
+{
+  r_8 vmin,vmax;
+  h2w_.GetMinMax(vmin,vmax);
+  double seuil=vmax/10000.;
+  if (seuil<1.e-6) seuil=1.e-6;
+  r_8 facnorm=((fgnorm1)?vmax:1.);
+  cout << " MultiDish::GetProjNoiseLevel Min,Max=" << vmin << " , " << vmax 
+       << " facnorm=" << facnorm << " seuil=" << seuil << endl;  
+  HProf hp(0,kmax_,nbin);
+  for(sa_size_t j=0; j<h2w_.NBinY(); j++) {
+    double y=h2w_.Y(j);
+    for(sa_size_t i=0; i<h2w_.NBinX(); i++) {
+      double x=h2w_.X(i);
+      double yw=h2w_(i,j);
+      if (yw<seuil) continue;
+      hp.Add(sqrt(x*x+y*y),facnorm/yw);
+    }
+  }
+  return hp;
+}
+
+HProf MultiDish::GetProjResponse(int nbin, bool fgnorm1)
+{
+  r_8 vmin,vmax;
+  h2w_.GetMinMax(vmin,vmax);
+  r_8 facnorm=((fgnorm1)?(1./vmax):1.);
+  cout << " MultiDish::GetProjResponse Min,Max=" << vmin << " , " << vmax 
+       << " facnorm=" << facnorm << endl;  
+  HProf hp(0,kmax_,nbin);
+  for(sa_size_t j=0; j<h2w_.NBinY(); j++) {
+    double y=h2w_.Y(j);
+    for(sa_size_t i=0; i<h2w_.NBinX(); i++) {
+      double x=h2w_.X(i);
+      hp.Add(sqrt(x*x+y*y),h2w_(i,j)*facnorm);
+    }
+  }
+  return hp;
+}
+
 
 Histo2D MultiDish::PosDist(int nx, int ny, double dmax)
@@ -312 +414 @@
       double kx0 = DeuxPI*(dishes_[i].X-dishes_[j].X)/lambda_;
       double ky0 = DeuxPI*(dishes_[i].Y-dishes_[j].Y)/lambda_;
+      double pgain=dishes_[i].Gain()*dishes_[j].Gain();
       rot.Do(kx0, ky0);
       //    if (kx0<0) kx0=-kx0;
@@ -321 +424 @@
           double y = jy*dy;
           if ((ix>0)&&(jy>0)) {
-            sumw += AddToHisto(kx0, ky0, x, y, rd(x,y), fgfh);
-            if (j!=i) sumw += AddToHisto(-kx0, -ky0, x, y, rd(x,y), fgfh);
+            sumw += AddToHisto(kx0, ky0, x, y, rd(x,y)*pgain, fgfh);
+            if (j!=i) sumw += AddToHisto(-kx0, -ky0, x, y, rd(x,y)*pgain, fgfh);
           }
           else {
             if ((ix==0)&&(jy==0)) {
-              sumw += h2w_.Add(kx0, ky0, rd(0.,0.), fgfh);
-              if (j!=i) sumw += h2w_.Add(-kx0, -ky0, rd(0.,0.), fgfh);
+              sumw += h2w_.Add(kx0, ky0, rd(0.,0.)*pgain, fgfh);
+              if (j!=i) sumw += h2w_.Add(-kx0, -ky0, rd(0.,0.)*pgain, fgfh);
             }
             else {
-              double w = rd(x,y);
+              double w = rd(x,y)*pgain;
               if (ix==0) {
                 sumw += h2w_.Add(kx0, ky0+y, w, fgfh);