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

/*  ------------------------ Projet BAORadio -------------------- 
    Classe  ForegroundCleaner
    R. Ansari , C. Magneville - Juin 2010 
---------------------------------------------------------------  */

#include "fgndsub.h"
#include "lobe.h"

#include "cubedef.h"
#include "matharr.h"

#include "ctimer.h"

/* --Methode-- */
ForegroundCleaner::ForegroundCleaner(Four2DResponse& arrep, Four2DResponse& tbeam, TArray< TF >& skycube)
  : arrep_(arrep) , tbeam_(tbeam), skycube_(skycube)
{
  double dxdeg = ThetaSizeDegre/(double)NTheta;
  double dydeg = PhiSizeDegre/(double)NPhi;
  dx_ = DegreeToRadian(dxdeg);
  dy_ = DegreeToRadian(dydeg);
  freq0_ = Freq0MHz;
  dfreq_ = FreqSizeMHz/(double)NFreq;
  cout << " ForegroundCleaner: " << " dx=" << dxdeg << " dy=" << dydeg << " degres ( dx_rad=" << dx_ << " dy_rad=" << dy_ << ")" 
       << " Freq0=" << freq0_ << " deltaFreq=" << dfreq_ << " MHz" << endl;
  skycube.Show();
}

/* --Methode-- */
void ForegroundCleaner::BeamCorrections()
{
  BeamEffect beam(arrep_);
  beam.Correct2RefLobe(tbeam_, skycube_, dx_, dy_, freq0_, dfreq_);
  cout << " ForegroundCleaner::BeamCorrections() done " << endl;
}

/* --Methode-- */
int ForegroundCleaner::CleanPointSources(double nsigmas)
{
  Timer tm("ForegroundCleaner::CleanPointSources");
  TArray< TF > sky2d(skycube_.SizeX(), skycube_.SizeY());
  for(sa_size_t ky=0; ky<sky2d.SizeY(); ky++) 
    for(sa_size_t kx=0; kx<sky2d.SizeX(); kx++) 
      sky2d(kx, ky) = skycube_(Range(kx), Range(ky), Range::all()).Sum();


  double mean, sigma;
  

  TArray< TF > amz(1,1,skycube_.SizeZ()), asz(1,1,skycube_.SizeZ());
  for(sa_size_t kz=0; kz<skycube_.SizeZ(); kz++)  {
    TArray< TF > slice = skycube_(Range::all() , Range::all(), Range(kz));
    MeanSigma(slice, mean, sigma);
    amz(0,0,kz)=mean;  asz(0,0,kz)=sigma; 
  }

  MeanSigma(sky2d, mean, sigma);
  cout << " ForegroundCleaner::CleanPointSources 2D Sky projection, mean=" << mean << " sigma=" << sigma << endl;  
  
  TF seuil = (TF)(mean+nsigmas*sigma);
  
  sa_size_t srccnt=0;
  for(sa_size_t ky=0; ky<skycube_.SizeY(); ky++) 
    for(sa_size_t kx=0; kx<skycube_.SizeX(); kx++) {
      if (sky2d(kx,ky)>seuil) {
        srccnt++; 
        skycube_(Range(kx), Range(ky), Range::all()) = amz;
      }
    }
  
  cout << " Cleaned NSrc= " << srccnt << " 2D source/pixels (TotNPix=" << sky2d.Size() 
       << ")-> " << 100.*srccnt/sky2d.Size() <<  "% with S>" << seuil << " NSigmas=" << nsigmas << endl;
  return (int)srccnt;
}


/* --Methode-- */
TArray< TF >  ForegroundCleaner::extractLSSCube(TArray< TF >& synctemp, TArray< TF >& specidx)
{
  Timer tm("ForegroundCleaner::extractLSSCube");
// Inputs : maplss, mapsyc, freq0, dfreq
// Outputs : synctemp, specidx  (reconstructed foreground temperature and spectral index
// Return_Array : foreground subtracted LSS signal 
{
  sa_size_t sz[5];   sz[0]=skycube_.SizeX();  sz[1]=skycube_.SizeY();
  synctemp.SetSize(2, sz); 
  specidx.SetSize(2, sz);
  TArray<r_4> omap; 
  omap.SetSize(skycube_, true);
  Vector vlnf(skycube_.SizeZ());
  int nprt = 0;
  // double freq0 : Frequence premier index en k (MHz)
  // double dfreq :   // largeur en frequence de chaque plan (Mhz)  
  for(sa_size_t i=0; i<skycube_.SizeX(); i++) 
    for(sa_size_t j=0; j<skycube_.SizeY(); j++)  {
      r_8 s1, sx, sx2, sy, sxy;
      s1=sx=sx2=sy=sxy=0.;
      for(sa_size_t k=0; k<skycube_.SizeZ(); k++)  {
        double lnf=log((double)k*dfreq_+freq0_);
        vlnf(k)=lnf;
        double ttot=(r_8)(skycube_(i,j,k));
        if (ttot < 1.e-5) continue;
        double lntt=log(ttot);
        s1+=1.;  sx+=lnf;  sx2+=(lnf*lnf);
        sy+=lntt;    sxy+=(lnf*lntt);
      }
      double beta = (sx*sxy-sx2*sy)/(sx*sx-s1*sx2);
      double alpha = (s1*sxy-sx*sy)/(s1*sx2-sx*sx);
      double T0 = exp(beta+alpha*vlnf(0));
      if ((i%16==0)&&(j%27==0)) 
        cout << "extractLSSCube[" << i << "," << j << "]: T0=" << T0 << " alpha=" << alpha 
             << " (mapsync=" << skycube_(i,j,0) << " ... " << skycube_(i,j,125) << ")" << endl; 
      synctemp(i,j) = T0;
      specidx(i,j) = alpha;
      for(sa_size_t k=0; k<skycube_.SizeZ(); k++) {
        r_4 fittedtemp = (r_4)(exp(beta+alpha*vlnf(k)));
        omap(i,j,k) = skycube_(i,j,k)-fittedtemp;
      }
    }
  return omap;
}

}