source: Sophya/trunk/Cosmo/RadioBeam/calcpk2.cc@ 4022

/*  ------------------------ Projet BAORadio -------------------- 
  Programme de calcul du spectre de puissance (3D) a partir d'un 
  cube de  delta T/T LSS, d'un cube  delta T/T LSS synchrotron
  ou radio-sources, apres ajustement / soustraction d'une loi de 
  puissance en frequence (l'axe Z du tableau doit etre en frequence)

    R. Ansari , C. Magneville - Juin 2010 

Usage: calcpk2 [-t -g -mxr val] InMapLSS convFacLSS InMapSync convFacSync InMapRadioSource convFacRsc OutPkFile 
               [PixNoiseLevel] [Diameter/Four2DRespTableFile] [TargetBeamArcmin] [NSigSrcThr]
---------------------------------------------------------------  */

#include "machdefs.h"
#include "sopnamsp.h"
#include <iostream>
#include <string>
#include <math.h>

#include <typeinfo>

#include "specpk.h"
#include "histats.h"
#include "vector3d.h"

#include "qhist.h" 
#include "lobe.h" 
#include "cubedef.h" 
#include "fgndsub.h" 
#include "radutil.h" 

#include "histinit.h"
#include "fftwserver.h"
#include "randr48.h"      

#include "ctimer.h"

typedef ThSDR48RandGen RandomGenerator ;

//-------------------------------------------------------------------------
//      ------------------ MAIN PROGRAM ------------------------------
//-------------------------------------------------------------------------
/* --Fonction-- */
int main(int narg, const char* arg[])
{
  if ( (narg<6)||((narg>1)&&(strcmp(arg[1],"-h")==0)) ) {
    cout << " Usage: [-t -g -mxr val -fmt T0,alpha] calcpk2 InMapLSS convFacLSS InMapFgnd convFacFgnd OutPkFile \n" 
         << "        [PixNoiseLevel] [D_Dish/Four2DRespTableFile CorBeamDiam] \n"
         << "        [NSigSrcThr] [P2/P1] [RecMapFile] " << endl;
    if ((narg>1)&&(strcmp(arg[1],"-h")==0)) {
      cout << "-t -g : Triangular / gaussian beam shape (def=gaussian) \n" 
           << "-fmt T0,alpha: fix mean slice temp. according to power law T0,alpha \n "
           << "-mxr val: Max beam correction factor (default=10.) \n "
           << "- InMapLSS: Input 3D LSS cube (PPF file name) \n "
           << "- convFacLSS: LSS cube conversion factor to mK (milliKelvin) \n"
           << "- InMapFgnd: Input 3D foreground cube (PPF file name) \n"
           << "- convFacFgnd: Foreground cube conversion factor to mK (milliKelvin) \n"
           << "- PixNoiseLevel: White noise level per pixel (mK) (default=0.) \n" 
           << "- D_Dish/Four2DRespTableFile: Dish diameter or 2D (u,v) plane response (PPF file name) \n"
           << "- CorBeamDiam: Beam correction target dish diameter  \n" 
           << "    These two parameters are used to correct for beam effect for a \n"
           << "    target beam (independent of frequency) defined by D/Lambda \n"
           << "    DoL = 100 --> beam ~ 35 arcmin (D=30m @ z~0.5 Lambda~30cm) \n" 
           << "    default : no beam correction applied \n " 
           << " - NSigSrcThr: Point source cleaning, Nb_Sigmas on stacked 2D temperature \n" 
           << "    default (0.) : no point source cleaning, use NSigSrcThr ~ 3..5 \n"
           << " - P2/P1: 2nd/first degree polynomial fit on ln(Temp) = f(ln(freq)) \n "
           << "    foreground subtraction. default is P2 \n" 
           << "- RecMapFile: output PPF file for reconstructed foreground template \n" 
           << "    (Temperature,SpectralIndex) and extracted LSS cube \n"
         << endl;
      return 1;
    }
    else cout << "   calcpk2 -h for detailed usage " << endl;
    return 2;
  }
  Timer tm("calcpk2");
  int rc = 0;
  try {
    bool fggaussian=true;  // true -> gaussian beam 
    double maxratio=10.;   // valeur max du rapport des lobes lors de la correction de lobe

    bool fgfix_mXYtemp=false;
    double T0_pl=1.; 
    double alpha_pl=0.;
    // decodage argument optionnel 
    bool fgoptarg=true;
    while (fgoptarg) {
      string fbo = arg[1];
      if (fbo=="-t")  { fggaussian=false; arg++; narg--; }
      else if (fbo=="-g")  { fggaussian=true; arg++; narg--; }
      else if (fbo=="-mxr")  { arg++; maxratio=atof(arg[1]); arg++; narg-=2; }
      else if (fbo=="-fmt")  
        { fgfix_mXYtemp=true; arg++; sscanf(arg[1],"%lg,%lg",&T0_pl,&alpha_pl); arg++; narg-=2; }
      else fgoptarg=false;
    }
    if (narg < 6) {
      cout << " calcpk2/error arguments , applobe -h for help " << endl;
      return 2;
    }

    string inppflss = arg[1];
    r_4 rfaclss = atof(arg[2]);
    string inppfsync = arg[3];
    r_4 rfacsync = atof(arg[4]);
    string outname = arg[5];

    double pixsignoise = 0.;
    bool fgaddnoise=false;
    if (narg>6) {
      pixsignoise=atof(arg[6]);
      if (pixsignoise>1.e-6)  fgaddnoise=true;
    }
    
    bool fgcorrbeam=true;

    bool fgresptbl=false;
    double DIAMETRE=100.;
    string resptblname;
    if (narg>7) {
      if (isdigit(*arg[7])) {
        fgresptbl=false;
        DIAMETRE=atof(arg[7]);
      }
      else { 
        resptblname=arg[7];
        fgresptbl=true;
      }
    }
    double tbeamDiam=0.;
    if (narg>8) {
      tbeamDiam=atof(arg[8]);
      if (tbeamDiam<1.)  fgcorrbeam=false;
    }
    bool fgclnsrc=true;
    double nsigsrc=5.;
    if (narg>9) {
      nsigsrc=atof(arg[9]);
      if (nsigsrc<1.e-6)  fgclnsrc=false;
    }
    bool fgpoly2=true;  // true -> soustraction polynome degre 2
    if ((narg>0)&&(strcmp(arg[10],"P1")==0)) fgpoly2=false;      
    bool fgsavemaps=false;
    string outmap_ppfname="extlss.ppf";
    if (narg>11) {
      outmap_ppfname=arg[11];
      fgsavemaps=true;
    }
    
    TArray<r_4> maplss, mapsync;
    const char * tits[2]={"LSS", "Sync/RadioSrc"};
    for(int ks=0; ks<2; ks++)  {
      string& ppfname=inppflss;
      r_4 rfac=rfaclss;
      TArray<r_4>* inmap=&maplss; 
      if (ks==1) {  ppfname=inppfsync;  rfac=rfacsync;  inmap=&mapsync; }
      cout << "calcpk2[" << ks+1 << "] : reading 3D map " << tits[ks] << " from file " << ppfname 
           << "  RenormFactor=" << rfac << endl;
      PInPersist pin(ppfname);
      pin >> (*inmap);
      (*inmap) *= rfac;
      double mean, sigma;
      MeanSigma(*inmap, mean, sigma);
      cout << " ...InMap sizes " << inmap->InfoString() << endl;
      inmap->Show(); 
      cout << " ... Mean=" << mean << " Sigma=" << sigma << endl;
    }

    bool smo; 
    if (!maplss.CompareSizes(mapsync,smo) ) {
      cout << " calcpk2/ERROR sizes " << endl;
      maplss.Show();  mapsync.Show(); 
      return 99;
    }
    
    TArray<r_4> skycube(mapsync);
    skycube += maplss;

    if (fgaddnoise) { 
      cout << " calcpk2: adding noise to skycube cube ... " << endl;
      BeamEffect::AddNoise(skycube, pixsignoise);
    }

    double mean, sigma;
    MeanSigma(skycube, mean, sigma);
    cout << " input sky cube : Mean=" << mean << " Sigma=" << sigma << endl;      
    tm.Split(" After input ");

    H21Conversions conv;
    conv.setFrequency(Freq0MHz);
    double lambda = conv.getLambda();
    Four2DResponse arep(2, DIAMETRE/lambda, DIAMETRE/lambda, lambda);
    Four2DResponse* arep_p=&arep;
    Four2DRespTable resptbl;
    if (fgresptbl) {
      cout << "calcpk2[3.a]: initializing Four2DRespTable from file" << resptblname << endl;
      resptbl.readFromPPF(resptblname);
      resptbl.renormalize(1.);
      arep_p=&resptbl;
    }
    else cout << " calcpk2[3.a]: Four2DResponse ( Diameter=" << DIAMETRE << " Lambda= " << lambda
              << " DoL=" << DIAMETRE/lambda << " ) " << endl;

    double DoL = tbeamDiam/lambda;
    double tbeamarcmin = RadianToDegree(1.22/DoL)*60.;
    int typcb = (fggaussian)?1:2;
    //    if (fgresptbl) typcb=22;
    Four2DResponse tbeam(typcb, DoL, DoL );

    ForegroundCleaner  cleaner(*arep_p, tbeam, skycube, maxratio); 
    if (fgfix_mXYtemp) {
      cout << "calcpk2[3.b] : calling cleaner.FixMeanXYTemp(T0=" << T0_pl << ",alpha=" << alpha_pl << ")" << endl;
      cleaner.FixMeanXYTemp(T0_pl,alpha_pl);
    }
    if (fgcorrbeam) {
      cout << "calcpk2[3.c] : calling cleaner.BeamCorrections() for target beam Diameter=" << tbeamDiam  
           << " D/Lambda=" << DoL  << "  -> arcmin " << tbeamarcmin << " TypDishResp=" << typcb << endl;
      cleaner.BeamCorrections();
    }
    cout << " calcpk2[3.d] : calling cleaner.CleanNegatives() ... " << endl;
    cleaner.CleanNegatives();
    if (fgclnsrc) {
      cout << "calcpk2[3.e] : calling cleaner.CleanPointSources() with threshold NSigma=" << nsigsrc  << endl;
      cleaner.CleanPointSources(nsigsrc);
    }

    cout << "calcpk2[4] : calling cleaner.extractLSSCube(...) " << endl;
    TArray<r_4> synctemp, specidx; 
    TArray<r_4> exlss; 
    if (fgpoly2) exlss = cleaner.extractLSSCubeP2(synctemp, specidx);
    else  exlss = cleaner.extractLSSCubeP1(synctemp, specidx);

    MeanSigma(exlss, mean, sigma);
    cout << " After cleaning/extractLSS: Mean=" << mean << " Sigma=" << sigma << endl;      
    tm.Split(" After CleanForeground");

    cout << "calcpk2[5] : computing 3D Fourier coefficients ... " << endl;
    FFTWServer ffts;
    TArray< complex<r_4> > four3d;
    ffts.FFTForward(exlss, four3d);
    tm.Split(" After FFTForward ");

    cout << "calcpk2[6] : computing power spectrum ... " << endl;
    RandomGenerator rg;
    Four3DPk pkc(four3d, rg);
    double dkxmpc = DeuxPI/(double)exlss.SizeX()/XCellSizeMpc;
    double dkympc = DeuxPI/(double)exlss.SizeY()/YCellSizeMpc;
    double dkzmpc = DeuxPI/(double)exlss.SizeZ()/ZCellSizeMpc;
    pkc.SetCellSize(dkxmpc, dkympc, dkzmpc);
    
    HProf hp = pkc.ComputePk(0.,HPk_NBin);

    tm.Split(" Done ComputePk ");    
    { 
    cout << "calcpk2[7.a] : writing profile P(k)  to  " << outname << endl;
    POutPersist po(outname);
    po << hp;
    }
    if (fgsavemaps) {
      cout << "calcpk2[7.b] : writing foreground maps and extracted LSS to  " << outmap_ppfname << endl;
      POutPersist pom(outmap_ppfname);
      pom << PPFNameTag("Tsync") << synctemp;
      pom << PPFNameTag("async") << specidx;
      pom << PPFNameTag("extlss") << exlss;
    }

  }  // End of try bloc 
  catch (PThrowable & exc) {  // catching SOPHYA exceptions
    cerr << " calcpk2.cc: Catched Exception (PThrowable)" << (string)typeid(exc).name() 
         << "\n...exc.Msg= " << exc.Msg() << endl;
    rc = 99;
  }
  catch (std::exception & e) {  // catching standard C++ exceptions
    cerr << " calcpk2.cc: Catched std::exception "  << " - what()= " << e.what() << endl;
    rc = 98;
  }
  catch (...) {  // catching other exceptions
    cerr << " calcpk2.cc: some other exception (...) was caught ! " << endl;
    rc = 97;
  }
  cout << " ==== End of calcpk2.cc program  Rc= " << rc << endl;
  return rc;    
}