source: Sophya/trunk/SophyaProg/PMixer/extractRS.cc@ 928

#include "machdefs.h"
#include <stdlib.h>
#include <stdio.h>
#include <fstream.h>
#include <iostream.h>
#include <math.h>

#include <string>
#include <vector>

#include "tmatrix.h"
#include "tvector.h"
#include "fioarr.h"
#include "vector3d.h"
#include "srandgen.h"
#include "timing.h"
#include "sambainit.h"
#include "pexceptions.h"
#include "datacards.h"
#include "fitsioserver.h"
#include "fitsspherehealpix.h"
#include "radspecvector.h"
#include "blackbody.h"
#include "nupower.h"

#include "squarefilt.h"
#include "trianglefilt.h"
#include "specrespvector.h"
#include "gaussfilt.h"
#include "fitstarray.h"


/*! \ingroup PMixer
    \file extractRS.cc
 
  \brief  \b PROGRAM \b extractRS  <BR>
  Program to create a map (SphereHEALPix) with point source
  from small map (matrix) and convolving their radiation spectrum
  RadSpectra
  with a given filter response SpectralResponse
*/

// -----------------------------------------------------------------
// ------------- Function declaration ------------------------------
int CheckCards(DataCards & dc, string & msg);
char * BuildFITSFileName(string const & fname);
SpectralResponse * getSpectralResponse(DataCards & dc);
template <class T>
void  MeanSig(NDataBlock<T> const & dbl, double& gmoy, double& gsig);
float ComputeFrequency(DataCards &,string&);

// -----------------------------------------------------------------

//  ----- Global (static) variables ------------
static bool rdmap = false;    // true -> Read map first 
static char mapPath[256];     // Path for input maps
static int  hp_nside = 32;    // HealPix NSide
static int  nskycomp = 0;     // Number of sky components
static int  debuglev = 0;     // Debug Level
static int  printlev = 0;     // Print Level
static POutPersist * so = NULL;  // Debug PPFOut file

// -------------------------------------------------------------------------
//                             main program 
// -------------------------------------------------------------------------
int main(int narg, char * arg[])
{
  if ((narg < 3) || ((narg > 1) && (strcmp(arg[1], "-h") == 0) )) {
    cout << "  Usage: extractRS parameterFile outputfitsnameformaps  [outppfname]" << endl;
    exit(0);
  }
  
  InitTim();
  
  string msg;
  int rc = 0;
  RadSpectra * es = NULL;
  SpectralResponse * sr = NULL;
  double moy, sig;
  
  DataCards dc;
  so = NULL;
  
  try {
    string dcard = arg[1];
    cout << " Decoding parameters from file " << dcard << endl;
    dc.ReadFile(dcard);
    
    rc = CheckCards(dc, msg);
    if (rc) { 
      cerr << " Error condition -> Rc= " << rc << endl;
      cerr << " Msg= " << msg << endl;
      return(rc);
    }
  }
  catch (PException exc) {
    msg = exc.Msg();
    cerr << " !!!! extractRS/ Readcard - Catched exception - Msg= " << exc.Msg() << endl;
    return(90);
  }   
  
  
  cout << " extractRS/Info : NComp = " <<  nskycomp << " SphereHEALPix_NSide= " << hp_nside << endl;
  cout << "  ... MapPath = " << (string)mapPath << "  DebugLev= " << debuglev 
       << "  PrintLev= " << printlev << endl;
  

  // We create an output persist file for writing debug objects
  if (debuglev > 0) so = new POutPersist("extrRSdbg.ppf");
  SphereHEALPix<float> SourceMap(hp_nside);
  SphereHEALPix<float> OutputMap(hp_nside);
  bool okout = true;
  
  try {
    
    //    FitsIoServer fios; // Our FITS IO Server
    char * flnm, buff[90];
    string key;
    string key2;
   
    double K = 1.;
    
    // Loop over sky frequencies components
    int maxOfSkyComp = nskycomp;
    int nbOfPix;
    int sk;
    Vector freq(maxOfSkyComp);
    TMatrix<float> spectrum;
    TMatrix<float> fdenu; //par freq et par pixel
    int nb_source;
    TMatrix<float> rareSource; //par freq et par numero de source
    
    for(sk = 0; sk<maxOfSkyComp; sk++) 
      {
        TMatrix<float> ings;    
        // Opening the file containing the map
        cout << "-------------------" << endl;
        cout << " Processing map's frequency component No " << sk+1 << endl;
        sprintf(buff, "%d", sk+1);
        key = (string)"RADSPECMAP" + buff;
        
        flnm = BuildFITSFileName(dc.SParam(key,0));
        cout << " Reading Input FITS map " << (string)flnm << endl;
        FITS_TArray<float> read_input(flnm); 
        ings = (TArray<float>)read_input;

        if (printlev > 2) 
          {
            double min = 1000000000;
            double max = -10000000000;
            for(int x=0;x<ings.NRows(); x++)
              {
                for(int y=0;y<ings.NCols(); y++)
                  {
                    if(min>ings(x,y)) min = ings(x,y);
                    if(max<ings(x,y)) max = ings(x,y);
                  }
              }
            cout << "min and max values for input map" << min << " " << max << endl;
          }
        // Computing the frequency according to the name of the file
        freq(sk) =  ComputeFrequency(dc,key);
        if (printlev > 2)  cout << "filling spectrum for freq = " << freq(sk) << endl;
        
        // Opening the file containing the point source
        sprintf(buff, "%d", sk+1);
        key2 = (string)"SOURCEPTMAP" + buff;
        flnm = BuildFITSFileName(dc.SParam(key2,0));
        const char* nameOfFile = flnm ;
        if (printlev > 2) cout << " Reading Input FITS map " << nameOfFile << endl;
        ifstream from(nameOfFile);
        char ch;
        string dum;
        from >> nb_source >> dum;
        if(sk==0) rareSource.ReSize(maxOfSkyComp,nb_source);
        for (int ii=0; ii< nb_source; ii++)
          {
            double value = 0;
            from >> value;
            rareSource(sk,ii) = value*1.e-26; //for Jansky->W
          }
        from.close();

        nbOfPix = ings.NRows()*ings.NCols();
        spectrum.ReSize(ings.NRows(),ings.NCols());
        if(sk==0) fdenu.ReSize(maxOfSkyComp,nbOfPix);
        int xy=0;
        for(int x=0;x<ings.NRows(); x++)
          {
            for(int y=0;y<ings.NCols(); y++)
              {
                spectrum(x,y) = ings(x,y)*1.e-26; //for Jansky->W
                fdenu(sk,xy) = spectrum(x,y);  
                xy += 1;
              }
          }
        if (printlev > 2) 
          {
            MeanSig(spectrum.DataBlock(), moy, sig );
            cout << " MeanSig for spectrum  map - Mean= " << moy << " Sigma= " << sig << endl;
          }
        
        
        K = dc.DParam(key, 1, 1.);
        if (debuglev > 4) {  // Writing the input map to the outppf
          FIO_TArray<float> fiog;
          fiog.Write(*so, key);
        }
        
        if(printlev>2) 
          {     
            sprintf(buff, "End of Proc. Comp. %d ", sk+1);
            cout << " --------------------------------------- \n" << endl;    
            PrtTim(buff);
          }
      }   // End of sky component loop 
    
    // ---------------------------------------------------------------------
    //      Integration in the detector band of the maps
    // ---------------------------------------------------------------------
    SpectralResponse* sr = NULL;
    sr = getSpectralResponse(dc);
  
    if (printlev > 2)  cout  << "SpectralResponse!" << *sr << endl;
    TVector<float> outCoeff(nbOfPix);

    for(int pix=0;pix<nbOfPix;pix++)
      {
        Vector smallFDeNu(maxOfSkyComp);
        for(sk = 0; sk<maxOfSkyComp; sk++) 
          {
            smallFDeNu(sk) =  fdenu(sk,pix);
          }
        
        RadSpectraVec radSV(freq,smallFDeNu);
        if (printlev > 10)  cout << "RadiationSpectra" << radSV << endl;

        outCoeff(pix) = radSV.filteredIntegratedFlux(*sr);
      }
    
    for(int i=0;i<OutputMap.NbPixels();i++)
      {
        float random = frand01();
        if(random == 1) random = frand01();
        int randomPix = random*OutputMap.NbPixels();
        int outpix = random*nbOfPix; 
        OutputMap(i) = outCoeff(outpix);
      }

    
    // ---------------------------------------------------------------------
    //      Integration in the detector band of the sources
    // ---------------------------------------------------------------------
    TVector<float> sourceCoeff(nb_source);
    for(int nsource=0;nsource<nb_source;nsource++)
      {
        Vector SourceFDeNu(maxOfSkyComp);
        for(sk = 0; sk<maxOfSkyComp; sk++) 
          {
            SourceFDeNu(sk) =  rareSource(sk,nsource);
          }
        RadSpectraVec radSV_Source(freq,SourceFDeNu);
        sourceCoeff(nsource) = radSV_Source.filteredIntegratedFlux(*sr);
      }

    for(int i=0;i<nb_source;i++)
      {
        float random = frand01();
        int outPix = random*OutputMap.NbPixels(); 
        SourceMap(outPix) = sourceCoeff(i); 
        OutputMap(outPix) = OutputMap(i)+ SourceMap(outPix);
      }
  }  
  catch (PException exc) {
    msg = exc.Msg();
    cerr << " !!!! extractRS - Catched PException  - Msg= " << exc.Msg() << endl;
    rc = 50;
  }   
  catch (PThrowable exc) {
    msg = exc.Msg();
    cerr << " !!!! extractRS - Catched PThrowable  - Msg= " << exc.Msg() << endl;
    rc = 50;
  }   
  catch (...) {
    cerr << " Une exception de type inconnue !!! " << endl;
    exit(99);
  }
  
  // Saving the output map in FITS format 
  if (okout) 
    {
      //       {
      //        FitsIoServer fios;
      //        fios.save(SourceMap, arg[3]);
      //        double moy,sig;
      //        MeanSig(SourceMap.DataBlock(),moy,sig);
      //        cout << " MeanSig for Source Map - Mean= " << moy << " Sigma= " << sig << endl;
      //        cout << "Source Map (SphereHEALPix<float>) written to FITS file " 
      //             << (string)(arg[3]) << endl;
      //       }
      {
         FITS_SphereHEALPix<float> fios2(OutputMap);
         fios2.Write(arg[2]);
         double moy,sig;
         MeanSig(OutputMap.DataBlock(),moy,sig);
         cout << " MeanSig for Output Map - Mean= " << moy << " Sigma= " << sig << endl;
         cout << "Output Map (SphereHEALPix<float>) written to FITS file " 
              << (string)(arg[2]) << endl;
      }
      if(printlev>2) PrtTim("End of WriteFITS ");
      // Saving the output map in PPF format 
      if (narg > 2) {
        POutPersist s(arg[3]); 
        FIO_SphereHEALPix<float> fiog(OutputMap);
        fiog.Write(s);
        cout << "Output Map (SphereHEALPix<float>) written to POutPersist file "  
             << (string)(arg[3]) << endl;
        if(printlev>2) PrtTim("End of WritePPF ");
      }
    }
  if (so) delete so;  //  Closing the debug ppf file 
  return(rc);
}



/* Nouvelle-Fonction */
int CheckCards(DataCards & dc, string & msg)
  //   Function to check datacards
{
  rdmap = false;
  mapPath[0] = '\0';
  nskycomp = 0;
  debuglev  = 0;
  printlev = 0;
  
  int rc = 0;
  string key, key2;
  
  //  Checking datacards
  if (dc.NbParam("EXT_RS") < 2)   {
     rc = 71; 
     msg = "Invalid parameters  - Check @EXT_RS card ";
     return(rc);
  }
  int kc;
  char buff[256];
  bool pb = false;
  string key3;
  int ncomp = 0;
  ncomp = dc.IParam("EXT_RS", 0);
  int nside = 32;
  nside = dc.IParam("EXT_RS", 1);
  
  for(kc=0; kc<ncomp; kc++) 
    {
      key = "RS_EXT_PATH";
      if (dc.NbParam(key) < 3)  strncpy(mapPath, dc.SParam(key, 0).c_str(), 255);
      mapPath[255] = '\0';
      sprintf(buff, "RADSPECMAP%d", kc+1);
      key = buff;
      
      if (dc.NbParam(key) < 1) {   
        msg = "Missing or invalid card : " + key;
        pb = true;  break;
      }
    }
  
  //  Checking detection filter specification 

  key = "GAUSSFILTER";
  key2 = "FILTERFITSFILE";
  if ( (dc.NbParam(key) < 5) && (dc.NbParam(key2) < 3)) { 
    msg = "Missing card or parameters : Check @GAUSSFILTER or @FILTERFITSFILE";
    rc = 73;  return(rc); 
  }
  
  if (pb)  {
    rc = 75;
    return(75);
  }
  
  
  // Initialiazing parameters
  rc = 0;
  msg = "OK";
  nskycomp = ncomp;
  hp_nside = nside;
  // Checking for PATH definition card
  key = "RS_EXT_PATH";
  if (dc.NbParam(key) < 3)  strncpy(mapPath, dc.SParam(key, 0).c_str(), 255);
  mapPath[255] = '\0';
  key = "DEBUGLEVEL";
  debuglev =  dc.IParam(key, 0, 0);
  key = "PRINTLEVEL";
  printlev =  dc.IParam(key, 0, 0);
  return(rc);
}

static char buff_flnm[1024];   // Mal protege !
/* Nouvelle-Fonction */
char* BuildFITSFileName(string const & fname)
{
if (mapPath[0] != '\0') sprintf(buff_flnm, "%s/%s", mapPath, fname.c_str());
else sprintf(buff_flnm, "%s", fname.c_str());
return(buff_flnm);
}


/* Nouvelle-Fonction */
template <class T>
void  MeanSig(NDataBlock<T> const & dbl, double& gmoy, double& gsig)

{
  gmoy=0.;
  gsig = 0.; 
  double valok;
  for(int k=0; k<(int)dbl.Size(); k++) {
    valok = dbl(k);
    gmoy += valok;  gsig += valok*valok; 
  }
  gmoy /= (double)dbl.Size();
  gsig = gsig/(double)dbl.Size() - gmoy*gmoy;
  if (gsig >= 0.) gsig = sqrt(gsig);
}

/* Nouvelle-Fonction */
SpectralResponse * getSpectralResponse(DataCards & dc)
{
  SpectralResponse * filt = NULL;

  string key = "FILTERFITSFILE";
  string key2 = "GAUSSFILTER";
  string ppfname = "filter";

  if (dc.HasKey(key) ) {      // Reading FITS filter file
    char ifnm[256];
    strncpy(ifnm, dc.SParam(key, 1).c_str(), 255);
    ifnm[255] = '\0';
    TMatrix<float> mtx(2,10);
    FITS_TArray<float> read_input(ifnm); 
    mtx = (TMatrix<float>)read_input;

    double numin = dc.DParam(key, 2, 1.);
    double numax = dc.DParam(key, 3, 9999.);
    Vector nu(mtx.NCols());
    Vector fnu(mtx.NCols());
    for(int k=0; k<mtx.NCols(); k++) {
      nu(k) = mtx(0, k);
      fnu(k) = mtx(1, k);
    }
  filt = new SpecRespVec(nu, fnu, numin, numax);
  ppfname = key;
  }
  else if (dc.HasKey(key2) ) {   // creating GaussianFilter 
    double nu0 = dc.DParam(key2, 0, 10.);
    double s = dc.DParam(key2, 1, 1.);
    double a = dc.DParam(key2, 2, 1.);
    double numin = dc.DParam(key2, 3, 0.1);
    double numax = dc.DParam(key2, 4, 9999);
    filt = new GaussianFilter(nu0, s, a, numin, numax);
    ppfname = key2;
    }
  if (filt == NULL) throw ParmError("datacard error ! No detection filter");
  cout << " Filter decoded - Created " << endl;
  cout << *filt << endl;

// for debug
  //  if (debuglev > 1)  SpectralResponse2Nt(*filt, *so, ppfname);
  return(filt);
}


float ComputeFrequency(DataCards & dc,string& key)
{
  ofstream filename("temp_filename");
  filename << dc.SParam(key,0) << '\n';
  filename.close();
  
  ifstream fromtempfile("temp_filename");
  char chtemp;
  string temp_str;
  int count = 0;
  do { 
    fromtempfile.get(chtemp);    
    if(count> 12 && count < 20)
      {
        temp_str+= chtemp; 
      }
    count+=1;
  } while(chtemp!='\n'); 
  fromtempfile.close();
  
  ofstream filename2("temp_filename");
  filename2 << temp_str << '\n';
  filename2.close();
  
  int frequency;
  ifstream fromtempfile2("temp_filename");
  fromtempfile2 >> frequency;
  fromtempfile2.close();
  return ((float)frequency/1000.);
}