source: Sophya/trunk/SophyaProg/PMixer/skymixer.cc@ 3891

#include "sopnamsp.h"
#include "pmixer.h"

/*!
 * \defgroup PMixer PMixer module
 * This module contains programs which:
 * <UL>
 * <LI> add several sky components, taking into account their
 * radiation spectra and convoluting them with a given filter
 * response : skymixer
 * <LI> create a map with point source : extractRS
 * <LI> generate sky components, radiation spectra and spectral 
 * response (small generator of maps) : tgsky and tgrsr
 * </UL> 
 * A detailed description may be found at:
 */
/*!
 * \ingroup PMixer
 * \file skymixer.cc
 *\brief   \b PROGRAM    \b skyMixer  <BR>
 *   add several sky components, taking into account their
 *radiation spectra and convoluting them with a given filter
 *response
 */

// -----------------------------------------------------------------
// ------------- Function declaration ------------------------------
int CheckCards(DataCards & dc, string & msg, FitsOutFile& fios);
char * BuildFITSFileName(string const & fname);
SpectralResponse * getSpectralResponse(DataCards & dc, FitsOutFile& fios);
RadSpectra * getEmissionSpectra(DataCards & dc, int nc, FitsOutFile& fios);
void SKM_MergeFITSKeywords(char * flnm);
void SKM_MergeFITSKeywords2(DVList & dvl);
void RadSpec2Nt(RadSpectra & rs, POutPersist & so, string name);
void SpectralResponse2Nt(SpectralResponse& sr, POutPersist & so, string name);
 
// to add different sky components and corresponding tools
//----------------------------------------------------------
template <class T> 
void addComponent(SpectralResponse&  sr, 
                                     PixelMap<T>& finalMap, 
                                     PixelMap<T>& mapToAdd, 
                                     RadSpectra& rs, double K=1.);
//
template <class T>
void addComponentBeta(SphereHEALPix<T>& finalMap, 
                      SphereHEALPix<T>& mapToAdd,SpectralResponse& sr, 
                      SphereHEALPix<T>& betaMap, double normFreq, double K);
//
template <class T>
void integratedMap(SpectralResponse&  sr,   
                   SphereHEALPix<T>& betaMap, double normFreq, SphereHEALPix<T>& intBetaMap);

//
template <class T>
void addComponentBeta(SphereHEALPix<T>& finalMap, 
                      SphereHEALPix<T>& mapToAdd,
                      SphereHEALPix<T>& intBetaMap, double K);
//
template <class T> 
void addDipole(SpectralResponse&  sr,  PixelMap<T>& finalMap, 
               double theta,double phi,double amp,double temp);
//
// -----------------------------------------------------------------

//  ----- 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
static DVList * dvl_fitskw = NULL; // Global DVList for all FITS Keywords 
bool checkpdmtype=false;
bool dvliston= false;
// --------- SkyMixer Version --------------
static double skm_version = 1.5;

// -------------------------------------------------------------------------
//                             main program 
// -------------------------------------------------------------------------
int main(int narg, char * arg[])
{
  if ((narg < 3) || ((narg > 1) && (strcmp(arg[1], "-h") == 0) )) {
    cout << "  Usage: skymixer parameterFile outputfitsname [outppfname]" << endl;
    exit(0);
  }
  
  InitTim();
  FitsOutFile fios(arg[2]);
     
  cout << " skymixer Version= " << skm_version << " SophyaVersion= " 
       << SophyaVersion() << endl;
  
  string msg;
  int rc = 0;
  RadSpectra * es = NULL;
  SpectralResponse * sr = NULL;
  double moy, sig;
  
  DataCards dc;
  so = NULL;
  // DVList for merging all FITS keywords and datacards
  dvl_fitskw = new DVList;
  
  try {
    string dcard = arg[1];
    if(printlev > 1) cout << " Decoding parameters from file " << dcard << endl;
    dc.ReadFile(dcard);
    
    rc = CheckCards(dc, msg,fios);
    if (rc) { 
      cerr << " Error condition -> Rc= " << rc << endl;
      cerr << " Msg= " << msg << endl;
      return(rc);
    }
  }
  catch (PException exc) {
    msg = exc.Msg();
    cerr << " !!!! skymixer/ Readcard - Catched exception - Msg= " << exc.Msg() << endl;
    return(90);
  }   
  
  
  cout << " skymix/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("skymixdbg.ppf");
  int countHead = 2;
  SphereHEALPix<float> outgs(hp_nside);
  try{
    if (rdmap) {  // Reading map from FITS file
      char ifnm[256];
      strncpy(ifnm, dc.SParam("READMAP", 0).c_str(), 255);
      ifnm[255] = '\0';
      cout << " Reading output HealPix map from FITS file " << (string)ifnm << endl;
      {
        FitsInFile fitsin(ifnm);
        fitsin >> outgs;
        //  Getting FITS keywords in primary header
        //      SKM_MergeFITSKeywords(ifnm); 
        //  Getting FITS keywords in data object header
        if(checkpdmtype)
          {
            int status = 1;
            string whatsup=fitsin.getStringKeyword(2,"PDMTYPE",status);
            if(!(status==0&&whatsup=="COMPMAP")) 
              {
                cerr << " !!!! skymixer - file" << ifnm << " is not a COMPMAP " << endl;
                throw ParmError("Error");
              }
            else cout << "read the comment COMPMAP from " << ifnm << endl;
          }
        if(dvliston) SKM_MergeFITSKeywords2(outgs.Info());
        else 
          {
            fios.appendInputHeader(fitsin, countHead);
          }
      }
      if(printlev>0)
        cout << " Output HealPIx Map read - NbPixels= " << 
          outgs.NbPixels() << endl;
      if (printlev > 0) {
        MeanSig(outgs.DataBlock(), moy, sig );
        cout << " MeanSig for outpout map - Mean= " << 
          moy << " Sigma= " << sig << endl;
      }
    }
    else {
      if(printlev>0)
        cout << " Output HealPix Map  created - NbPixels= " << 
          outgs.NbPixels() << endl;
      outgs.SetPixels(0.);
    }
    
    // Decoding detection pass-band filter 
    sr = getSpectralResponse(dc,fios);
    PrtTim(" After FilterCreation ");
    
    char * flnm, buff[90];
    string key;
    
    double K = 1.;
    double freqOfMap = 1.;
    // Loop over sky component 
    int sk;
    for(sk = 0; sk<nskycomp; sk++) {
      cout << "------------------------------------" << endl;
      cout << " Processing sky component No " << sk+1 << endl;
      

      sprintf(buff, "%d", sk+1);
      key = (string)"DIPOLE" + buff;
      // check for an eventual dipole
      if(dc.HasKey(key)) 
        {
          if (es) { delete es; es = NULL; }
          double temp = -10.;
          double theta=  dc.DParam(key,1,1.);
          double phi  =  dc.DParam(key,2,1.);
          double amp  =  dc.DParam(key,3,1.);
          if(dc.NbParam(key)>3) 
            {
              temp =  dc.DParam(key,4,1.);
            }
          cout << " creating dipole " << temp << " " << amp << " " << phi << " " << theta << " " << endl;
          addDipole(*sr, outgs,theta,phi,amp,temp);
        }
      else
        {
          sprintf(buff, "%d", sk+1);
          key = (string)"MAPFITSFILE" + buff;
          flnm = BuildFITSFileName(dc.SParam(key, 0));
          
          K = dc.DParam(key, 1, 1.);
          
          cout << " Reading Input FITS map " << (string)flnm << endl;
          SphereHEALPix<float> ings(hp_nside);
          {
            FitsInFile fiosIn(flnm);
            fiosIn >> ings;
            //  Getting FITS keywords in primary header
            //  SKM_MergeFITSKeywords(flnm); 
            //  Getting FITS keywords in data object header
            //if(fiosIn.hasKeyword(string("PDMTYPE"), 1)) cout << "read the comment COMPMAP from " << flnm << endl;
            if(checkpdmtype)
              {
                int status = 1;
                string ch = "PDMTYPE";
                string whatsup=fiosIn.getStringKeyword(2,ch,status);
                if(!(status==0&&whatsup=="COMPMAP")) 
                  {
                    cerr << " !!!! skymixer - file" << flnm << " is not a COMPMAP " << endl;
                    throw ParmError("Error");
                  }
                else cout << "read the comment COMPMAP from " << flnm << endl;
              }
            if(dvliston) SKM_MergeFITSKeywords2(ings.Info());
            else 
              {
                fios.appendInputHeader(fiosIn, countHead);
              }
          }
          if (debuglev > 4) {  // Writing the input map to the outppf
            FIO_SphereHEALPix<float> fiog(ings);
            fiog.Write(*so, key);
          }
          if (printlev > 2) {
            MeanSig(ings.DataBlock(), moy, sig );
            cout << " MeanSig for input map - Mean= " << moy << " Sigma= " << sig << endl;
          }
          bool mapDependentOfFreq = false;
          key = (string)"BETAFITSFILE"+ buff;
          if(dc.HasKey(key)) 
            {
              mapDependentOfFreq = true;
            }
          
          // getting Emission spectra  
          if(!mapDependentOfFreq)
            {
              if (es) { delete es; es = NULL; }
              es = getEmissionSpectra(dc, sk,fios);
              addComponent(*sr, outgs, ings, *es, K);
            }
          else
            {
              key = (string)"BETAFITSFILE"+ buff;
              //SphereHEALPix<float> betaMap;
              flnm = BuildFITSFileName(dc.SParam(key, 0));
              double normFreq = dc.DParam(key, 1, 1.);
              if (printlev > 4) cout << "....BetaFits... normalization Freq = " << normFreq << endl;
              int nSideForInt = dc.DParam(key, 2, 1.);
              if (printlev > 4) cout << "....BetaFits... NSide for Integration map = " << nSideForInt << endl;
              cout << "....BetaFits...  Reading Beta FITS map " << (string)flnm << endl;
              SphereHEALPix<float> betaMap(hp_nside);
              {
                FitsInFile fiosBM(flnm);
                fiosBM >> betaMap;
                //  Getting FITS keywords in primary header
                // SKM_MergeFITSKeywords(flnm); 
                //  Getting FITS keywords in data object header
                //if(fiosBM.hasKeyword(string("COMPMAP"), 1)) cout << "read the comment COMPMAP from " << flnm << endl;
                if(checkpdmtype)
                  {
                    int status = 1;
                    string whatsup=fiosBM.getStringKeyword(2,"PDMTYPE",status);
                    if(!(status==0&&whatsup=="INDEXMAP")) 
                      {
                        cerr << " !!!! skymixer - file" << flnm << " is not an INDEXMAP " << endl;
                        throw ParmError("Error");
                      }
                    else cout << "read the comment INDEXMAP from " << flnm << endl;
                  }
                if(dvliston) SKM_MergeFITSKeywords2(betaMap.Info()); 
                else 
                  {
                    fios.appendInputHeader(fiosBM, countHead);
                  }
              }
              if (printlev > 2) {
                MeanSig(betaMap.DataBlock(), moy, sig );
                cout << " MeanSig for Beta map - Mean= " << moy << " Sigma= " << sig << endl;
              }
              if (debuglev > 4) {  // Writing the input map to the outppf
                FIO_SphereHEALPix<float> fiogs(betaMap);
                fiogs.Write(*so, key);
              }
              
              if(nSideForInt<0) nSideForInt = sqrt((double)betaMap.NbPixels()/12);
              bool bydefault = true;
              if(!bydefault)
                addComponentBeta(outgs,ings,*sr,betaMap,normFreq, K);
              else
                {
                  // integrate the betamap over the SpectralResponse
                  SphereHEALPix<float> intBetaMap(nSideForInt);
                  integratedMap(*sr, betaMap, normFreq, intBetaMap);
                  if (debuglev > 4) {  // Writing the input map to the outppf
                    FIO_SphereHEALPix<float> fiogs2(intBetaMap);
                    fiogs2.Write(*so, "INTBETAMAP");
                  }
                  
                  betaMap.Resize(8);
                  if (printlev > 4)
                    {
                      MeanSig(intBetaMap.DataBlock(), moy, sig );
                      cout << "....BetaFits... MeanSig for intBetaMap - Mean= " 
                           << moy << " Sigma= " << sig << endl;
                    }
                  // add the integrated beta map 
                  addComponentBeta(outgs,ings,intBetaMap, K);
                }
            }
          
          MeanSig(outgs.DataBlock(), moy, sig );
          cout << " MeanSig for Sum map - Mean= " << moy << " Sigma= " << sig << endl;
          cout << "-------------------------------------------------" << endl;
          
          sprintf(buff, "End of Processing Component %d ", sk+1);
          PrtTim(buff);
        } 
    }
  }
  catch(PException exc) 
    {
      cout << "catched PException" << endl;
      msg = exc.Msg();
      cerr << " !!!! skymixer - Catched exception - Msg= " << exc.Msg() << endl;
      rc = 50;
      return(50);
    } 
  
 try {
 // Saving the output map in FITS format 
   cout << "Output Map (SphereHEALPix<float>) written to FITS file " 
        << (string)(arg[2]) << endl;
   {
     //     FitsOutFile fios(arg[2]);
     fios.firstImageOnPrimaryHeader(false); // Use secondary header 
     DVList& dvl = (*dvl_fitskw);
     dvl["PDMTYPE"] = "COMPMAP";
     dvl.SetComment("PDMTYPE", "Planck Data Model Type"); 
     dvl["SOPHYVER"] =  SophyaVersion();
     dvl.SetComment("SOPHYVER", "Sophya Version number"); 
     dvl["SKYMVER"] =  skm_version;
     dvl.SetComment("SKYMVER", "skymixer Version number"); 
     //     fios.DVListIntoPrimaryHeader(dvl);  
     dvl["CREATOR"] = "SkyMixer";
     outgs.Info() = dvl;
     fios << outgs ;
   }
   PrtTim("End of WriteFITS ");
   // Saving the output map in PPF format 
   if (narg > 3) {
     POutPersist s(arg[3]); 
     FIO_SphereHEALPix<float> fiog(&outgs) ;
     fiog.Write(s);
     cout << "Output Map (SphereHEALPix<float>) written to POutPersist file "  
          << (string)(arg[3]) << endl;
     PrtTim("End of WritePPF ");
   }
 }
 catch(PException exc) 
   {
     cout << "catched PException (2)" << endl;
     msg = exc.Msg();
     cerr << " !!!! skymixer(2) - Catched exception - Msg= " << exc.Msg() << endl;
     rc = 55;
     return(55);
   } 

 if (so) delete so;  //  Closing the debug ppf file 
 return(rc);
}

/* Nouvelle-Fonction */
int CheckCards(DataCards & dc, string & msg, FitsOutFile& fios)
//   Function to check datacards
{
rdmap = false;
mapPath[0] = '\0';
hp_nside = 32;
nskycomp = 0;
debuglev  = 0;
printlev = 0;

int rc = 0;
string key, key2,key3;

 DVList & dvl = *dvl_fitskw; 
  //  Cheking datacards
  if (dc.NbParam("SKYMIX") < 2)   {
     rc = 71; 
     msg = "Invalid parameters  - Check @SKYMIX card ";
     return(rc);
  }
  if(dvliston)
    dvl.SetS("SKYMIX", dc.GetParams("SKYMIX"));
  else
    {
      fios.insertCommentLineOnHeader("--------------------");
      fios.insertCommentLineOnHeader("SkyMixer DataCard");
      fios.insertCommentLineOnHeader("--------------------");
 
       string key = "SKYMIX"+ (string)dc.GetParams("SKYMIX");
       fios.insertCommentLineOnHeader(key);
    }

  key = "CHECKPDMTYPE";
  if (dc.HasKey(key)) {
    checkpdmtype=true;
  }

  
  key = "READMAP";
  if (dc.HasKey(key)) {
    if (dc.NbParam(key) < 1) {
     rc = 72; 
     msg = "Invalid parameters  - Check @READMAP card ";
     return(rc);
    }
    else rdmap = true;
    
    if(dvliston) dvl.SetS(key, dc.GetParams(key));
    else 
      {
        string keynew = key + (string)dc.GetParams(key);
        fios.insertCommentLineOnHeader(keynew);
      }
  }

//  Checking detection filter specification 
  key = "GAUSSFILTER";
  key2 = "FILTERFITSFILE";
  key3 = "DIPOLE";
  if ( (dc.NbParam(key) < 5) && (dc.NbParam(key2) < 3) && (dc.NbParam(key3) < 3)) { 
    msg = "Missing card or parameters : Check @GAUSSFILTER or @FILTERFITSFILE or @DIPOLE";
    rc = 73;  return(rc); 
    }

  if (dc.HasKey(key))   
    { 
      if(dvliston) dvl.SetS("GAUSFILT", dc.GetParams(key));
      else 
        {
          string keynew = key + (string)dc.GetParams(key);
          fios.insertCommentLineOnHeader(keynew);
        }
    }
  if (dc.HasKey(key2))
    {
      if(dvliston) dvl.SetS("FILTFILE", dc.GetParams(key2));
      else
        {
          string keynew = key2 + (string)dc.GetParams(key2);
          fios.insertCommentLineOnHeader(keynew);
        }
    }
  if (dc.HasKey(key3))  
    {
      if(dvliston) dvl.SetS("DIPOLE", dc.GetParams(key3));
      else
        {
          string keynew = key3 + (string)dc.GetParams(key3);
          fios.insertCommentLineOnHeader(keynew);
        }
    }
  // Decoding number of component and pixelisation parameter
  int mg = 32;
  int ncomp = 0;
  ncomp = dc.IParam("SKYMIX", 0, 0);
  mg = dc.IParam("SKYMIX", 1, 32);
  if (ncomp < 1)  {
    msg = "Invalid parameters  - Check datacards @SKYMIX ";
    rc = 74;
    return(rc);
  }

  // Checking detection filter specification
  //  Checking input FITS file specifications 
  int kc;
  char buff[256];
  bool pb = false;
  string key4;
  string key5;
  string key6;
  for(kc=0; kc<ncomp; kc++) {
    sprintf(buff, "MAPFITSFILE%d", kc+1);
    key = buff;
    sprintf(buff, "DIPOLE%d", kc+1);
    key3 = buff;
    if (dc.NbParam(key) < 1 && (dc.NbParam(key3)<1)) 
      {   
        msg = "Invalid card or Missing card" + key + " or " + key3 ;
        pb = true;  break;
      }
        
    sprintf(buff, "MAPFIL%d", kc+1);
    if (dc.HasKey(key))  
      {
        if(dvliston) dvl.SetS(buff, dc.GetParams(key));
        else 
          {
            string keynew = buff + (string)dc.GetParams(key);
            fios.insertCommentLineOnHeader(keynew);
          }
      } 
    sprintf(buff, "DIPOLE%d", kc+1);
    if (dc.HasKey(key3))  
      {
        if(dvliston) dvl.SetS(buff, dc.GetParams(key3));
        else
          {
            string keynew = buff + (string)dc.GetParams(key3);
            fios.insertCommentLineOnHeader(keynew);
          }
      }
    sprintf(buff, "SPECTRAFITSFILE%d", kc+1);
    key = buff;
    sprintf(buff, "BLACKBODY%d", kc+1);
    key2 = buff;
    sprintf(buff, "POWERLAWSPECTRA%d", kc+1);
    key3 = buff;
    sprintf(buff, "BETAFITSFILE%d", kc+1);
    key4 = buff;
    sprintf(buff, "DIPOLE%d", kc+1);
    key5 = buff;
    sprintf(buff, "DERIVBB%d", kc+1);
    key6 = buff;
    if ( (dc.NbParam(key) < 3) && (dc.NbParam(key2) < 1) && (dc.NbParam(key3) < 6) && (dc.NbParam(key4)<2)
         &&  (dc.NbParam(key6)<1) &&  (dc.NbParam(key5)<1))  {   
      msg = "Missing card or invalid parameters : " + key + " " + key2 + " " + key3 + " " + key4+ " " + key5;
      pb = true;  break;
    }

    sprintf(buff, "SPECTF%d", kc+1);
    if (dc.HasKey(key)) {   
      if(dvliston) dvl.SetS(buff, dc.GetParams(key));
      else { 
        string keynew = buff + (string)dc.GetParams(key);
        fios.insertCommentLineOnHeader(keynew);}
    }
    sprintf(buff, "BLACKB%d", kc+1);
    if (dc.HasKey(key2)) {
      if(dvliston) dvl.SetS(buff, dc.GetParams(key2));
      else { 
        string keynew = buff + (string)dc.GetParams(key2);
        fios.insertCommentLineOnHeader(keynew);}
    }

    sprintf(buff, "PLAWSP%d", kc+1);
    if (dc.HasKey(key3)) {
      if(dvliston) dvl.SetS(buff, dc.GetParams(key3));
      else { 
        string keynew = buff + (string)dc.GetParams(key3);
        fios.insertCommentLineOnHeader(keynew);}
    }

    sprintf(buff, "BETAFI%d", kc+1);
    if (dc.HasKey(key4))  {
      if(dvliston) dvl.SetS(buff, dc.GetParams(key4));
      else { 
        string keynew = buff + (string)dc.GetParams(key4);
        fios.insertCommentLineOnHeader(keynew);}
    }

    sprintf(buff, "DIPOLE%d", kc+1);
    if (dc.HasKey(key5))  {
      if(dvliston) dvl.SetS(buff, dc.GetParams(key5));
      else { 
        string keynew = buff + (string)dc.GetParams(key5);
        fios.insertCommentLineOnHeader(keynew);}
    }

    sprintf(buff, "DERIBB%d", kc+1);
    if (dc.HasKey(key6))  {
      if(dvliston) dvl.SetS(buff, dc.GetParams(key6));
      else { 
        string keynew = buff + (string)dc.GetParams(key6);
        fios.insertCommentLineOnHeader(keynew);}
    }
    
  }

  if (pb)  {
    rc = 75;
    return(75);
  }


// Initialiazing parameters
  rc = 0;
  msg = "OK";
  nskycomp = ncomp;
  hp_nside = mg;  


// Checking for PATH definition card
  key = "MAPPATH";
  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 */
SpectralResponse * getSpectralResponse(DataCards & dc, FitsOutFile& fios)
{
  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';
    Matrix mtx;
    FitsInFile fiis(ifnm);
    fiis.firstImageOnPrimaryHeader(false); // Use secondary header HDU=2
    fiis >> mtx ; 
    //  Getting FITS keywords in primary header
    //    SKM_MergeFITSKeywords(ifnm); 
    //  Getting FITS keywords in data object header
    if(checkpdmtype)
      {
        int status = 1;
        string whatsup=fiis.getStringKeyword(2,"PDMTYPE",status);
        if(!(status==0&&whatsup=="COMPMAP")) 
          {
            cerr << " !!!! skymixer - file" << ifnm << " is not a COMPMAP " << endl;
            throw ParmError("Error");
          }
        else cout << "read the comment COMPMAP from " << ifnm << endl;
      }
    if(dvliston) SKM_MergeFITSKeywords2(mtx.Info()); 
    else 
      {
        fios.appendInputHeader(fiis, 2);
      }
    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");
  if(printlev>0)
    {
      cout << endl;
      cout << " Filter decoded - Created " << endl;
      cout << *filt << endl;
    }
  // for debug
  if (debuglev > 1)  SpectralResponse2Nt(*filt, *so, ppfname);
  return(filt);
}

/* Nouvelle-Fonction */
RadSpectra * getEmissionSpectra(DataCards & dc, int nc, FitsOutFile& fios)
{
  char numb[16];
  sprintf(numb, "%d", nc+1);

  string key = (string)"SPECTRAFITSFILE" + numb;
  string key2 = (string)"BLACKBODY" + numb;
  string key5 = (string)"DERIVBB" + numb;
  string key3 = (string)"POWERLAWSPECTRA" + numb;
  string ppfname = "espectra";

  RadSpectra * rs = NULL;
  if (dc.HasKey(key) ) {    // Reading emission spectra from file
    char * ifnm = BuildFITSFileName(dc.SParam(key, 0));
    cout << " Reading Input FITS spectra file " << (string)ifnm << endl;
    Matrix mtx;
    FitsInFile fiis(ifnm);
    fiis.firstImageOnPrimaryHeader(false); // Use secondary header HDU=2
    fiis >> mtx ; 
    //  Getting FITS keywords in primary header
    //    SKM_MergeFITSKeywords(ifnm); 
    //  Getting FITS keywords in data object header
    //if(fiis.hasKeyword(string("PDMTYPE"), 1)) cout << "read the comment COMPMAP from " << ifnm << endl;
    if(checkpdmtype)
      {
        int status = 1;
        string whatsup=fiis.getStringKeyword(2,"PDMTYPE",status);
        if(!(status==0&&whatsup=="FGRSPEC")) 
          {
            cerr << " !!!! skymixer - file" << ifnm << " is not a  FGRSPEC" << endl;
            throw ParmError("Error");
          }
        else cout << "read the comment FGRSPEC from " << ifnm << endl;
      }
    if(dvliston) SKM_MergeFITSKeywords2(mtx.Info()); 
    else 
      {
        fios.appendInputHeader(fiis, 2);
      }
    double numin = dc.DParam(key, 2, 1.);
    double numax = dc.DParam(key, 3, 9999.);
    Vector nu(mtx.NCols());
    Vector tnu(mtx.NCols());
    for(int k=0; k<mtx.NCols(); k++) {
      nu(k) = mtx(0, k);
      tnu(k) = mtx(1, k);
    }
  rs = new RadSpectraVec(nu, tnu, numin, numax);
  ppfname = key;
  }
  else if (dc.HasKey(key2) ) { // Creating BlackBody emission spectra
    rs = new BlackBody(dc.DParam(key2, 0, 2.726));
    ppfname = key2;
  }
  else if (dc.HasKey(key5) ) { // Creating Dipole
    rs = new DerivBlackBody(dc.DParam(key5, 0, 3.E-3));
    ppfname = key5;
  }
  else if (dc.HasKey(key3) ) { // Creating PowerLaw emission spectra
    double a = dc.DParam(key3, 0, 1.);
    double nu0 = dc.DParam(key3, 1, 100.);
    double dnu = dc.DParam(key3, 2, 10.);
    double b = dc.DParam(key3, 3, 0.);
    double numin = dc.DParam(key3, 4, 0.1);
    double numax = dc.DParam(key3, 5, 9999);
    rs = new PowerLawSpectra(a, b, nu0, dnu, numin, numax);
    ppfname = key3;
  }
  if (rs == NULL) throw ParmError("datacard error ! missing Emission spectra");
  cout << " Emission spectra decoded - Created (" << ppfname << ")" << endl;
  cout << *rs << endl;
// for debug
  if (debuglev > 2)  RadSpec2Nt(*rs, *so, ppfname);
  return(rs);
}




template <class T> 
void addDipole(SpectralResponse&  sr,  PixelMap<T>& finalMap, 
               double theta,double phi,double amp,double temp)
{
  DerivBlackBody dbb;
  if(temp>0) dbb.setTemperature(temp);
  double coeff = dbb.filteredIntegratedFlux(sr) * amp;
  UnitVector vd(theta,phi);
  UnitVector vc(theta,phi);
  
  for(int i=0; i<finalMap.NbPixels(); i++)
    {
      double thetar,phir;
      finalMap.PixThetaPhi(i,thetar,phir);      
      vc.SetThetaPhi(thetar,  phir);
      finalMap(i) += vd.Psc(vc)*coeff;
    }
  if (debuglev > 4) {  // Writing the input map to the outppf
    SphereHEALPix<float> ings((int)sqrt((double)finalMap.NbPixels()/12));
    for(int i=0; i<finalMap.NbPixels(); i++)
      {
        double thetar,phir;
        finalMap.PixThetaPhi(i,thetar,phir);      
        vc.SetThetaPhi(thetar,  phir);
        ings(i) = vd.Psc(vc);
      }
    FIO_SphereHEALPix<float> fiog(ings);
    fiog.Write(*so, "dipole");
    cout << "Debug the dipole map....saved in debug file !" << endl;
  }
}
/* Nouvelle-Fonction */
template <class T> 
void addComponent(SpectralResponse&  sr, PixelMap<T>& finalMap, 
                  PixelMap<T>& mapToAdd, RadSpectra& rs, double K)
{
  // finalMap = finalMap + coeff* mapToAdd
  // coeff    =  convolution of sr and rs 
  // compute the coefficient corresponding to mapToAdd
  if (finalMap.NbPixels() != mapToAdd.NbPixels())   
    throw SzMismatchError("addComponent()/Error: Unequal number of Input/Output map pixels");
  double coeff = rs.filteredIntegratedFlux(sr) * K;
  if (printlev > 1) 
    cout << " addComponent - Coeff= " << coeff << " (K= " << K << ")" << endl; 
  for(int i=0; i<finalMap.NbPixels(); i++)
    {
      finalMap(i) += coeff * mapToAdd(i);
    }
}
/* Nouvelle-Fonction */
template <class T> 
void addComponentBeta(SphereHEALPix<T>& finalMap, 
                      SphereHEALPix<T>& mapToAdd,SpectralResponse& sr, 
                      SphereHEALPix<T>& betaMap, double normFreq, double K)
{
  // finalMap = finalMap + coeff* mapToAdd
  // coeff    =  convolution of sr and rs 
  // compute the coefficient corresponding to mapToAdd
  // betaMap is the map of (beta(theta,phi)) 

  int nbpix = finalMap.NbPixels();
  if (nbpix != mapToAdd.NbPixels())   
    throw SzMismatchError("addComponentBeta()/Error: Unequal number of Input/Output map pixels");
  if (printlev > 1) 
    {
      cout << "addComponentBeta - Coeff= "  << K  << endl; 
      cout << "nb pixels: " << finalMap.NbPixels() << endl;
    }
  SphereHEALPix<T> bigBetaMap((int)sqrt((double)nbpix/12));
  if(nbpix != betaMap.NbPixels())
    {
       Sph2Sph(betaMap,bigBetaMap);
     }
  for(int i=0; i<finalMap.NbPixels(); i++)
    {
      // coeff = integration of (nu/normFreq)^(-beta(theta,phi)) in each pixels
      RadSpectra* rs =  new PowerLawSpectra
        (1.,-bigBetaMap(i), 0., normFreq, 0.01, 800.);
      double coeff = rs->filteredIntegratedFlux(sr);
      finalMap(i) += coeff*K*mapToAdd(i);
    }
}

template <class T>
void integratedMap(SpectralResponse&  sr,   
                   SphereHEALPix<T>& betaMap, 
                   double normFreq, 
                   SphereHEALPix<T>& intBetaMap)
{
  PowerLawSpectra rs(1.,-2., 0., normFreq);
  
  if(betaMap.NbPixels()!=intBetaMap.NbPixels()) 
    {
      Sph2Sph(betaMap,intBetaMap);
      for(int i=0; i<intBetaMap.NbPixels(); i++)
        {
          rs.setExp(-intBetaMap(i));
          double coeff = rs.filteredIntegratedFlux(sr);
          intBetaMap(i) = coeff;
        }
    }
  else
    {      
      for(int i=0; i<intBetaMap.NbPixels(); i++)
        {
          rs.setExp(-betaMap(i));
          double coeff = rs.filteredIntegratedFlux(sr);
          intBetaMap(i) = coeff;
        }
    }
}

template <class T>
void addComponentBeta(SphereHEALPix<T>& finalMap, 
                      SphereHEALPix<T>& mapToAdd,SphereHEALPix<T>& intBetaMap, double K)
{
  // finalMap = finalMap + coeff* mapToAdd
  // coeff    =  convolution of sr and rs 
  // compute the coefficient corresponding to mapToAdd
  // integBetaMap is the map of the integration (nu/normFreq)^(-beta(theta,phi)) over
  // the spectralResponse
  // different from addComponentBeta(PixelMap<T>& finalMap, 
  //    PixelMap<T>& mapToAdd,SpectralResponse& sr, PixelMap<T>& betaMap, double normFreq, double K)
  //   since it permits to use a intBetaMap with a different number of pixels than
  //   the other maps

  int nbpix = finalMap.NbPixels();
  if (nbpix != mapToAdd.NbPixels())   
    throw SzMismatchError("addComponentBeta(PixelMap<T>&,PixelMap<T>&,PixelMap<T>&,double)/Error: Unequal number of Input/Output map pixels");
  double coeff = K;

  if(nbpix != intBetaMap.NbPixels())
    {
      for(int i=0; i<finalMap.NbPixels();i++)
        {
          double teta,phi,val;
          finalMap.PixThetaPhi(i, teta, phi);
          int pixel = intBetaMap.PixIndexSph(teta,phi);
          val = intBetaMap.PixVal(pixel);
          finalMap(i) += coeff*mapToAdd(i)*val;
        }
    }
  else
    {
      for(int i=0; i<finalMap.NbPixels();i++)
        {
          finalMap(i) += coeff*mapToAdd(i)*intBetaMap(i);
        }
    }
  if (printlev > 1) 
    {
      cout << "addComponentBeta(SG<T>,SG<T>,SG<T>,double) - Coeff= "  << K  << endl;     
    }
}





/* Nouvelle-Fonction */
void SKM_MergeFITSKeywords(char * flnm)
{
  DVList dvl;
  FitsFile::FitsExtensionType typeOfExtension;
  int naxis;
  vector<int> naxisn;
  FitsFile::FitsDataType dataType;
  FitsInFile::GetBlockType(flnm, 1, typeOfExtension, naxis, naxisn, dataType, dvl);
  // Cleaning the keywords 
#define SZexlst 21
  const char *exlst[SZexlst]=
  {"SIMPLE","BITPIX" ,"NAXIS" ,"NAXIS#" ,"PCOUNT","GCOUNT",
     "EXTEND","ORIGIN" ,"DATE*" ,"TFIELDS","TTYPE#","TFORM#",
     "TUNIT#","EXTNAME","CTYPE#","CRVAL#" ,"CRPIX#","CDELT#",
     "XTENSION","INSTRUME","TELESCOP"};
  char  kwex[32];
  int i,l;
  for (i=0; i<SZexlst; i++) {
    strncpy(kwex, exlst[i], 32);
    l = strlen(kwex)-1;
    if ((kwex[l] != '*') && (kwex[l] != '#')) {
      dvl.DeleteKey(kwex);
    } 
    else {
      bool fgd = (kwex[l] == '#') ? true : false;
      list<string> lstsup;
      kwex[l] = '\0';
      DVList::ValList::const_iterator it;
      for (it = dvl.Begin(); it != dvl.End(); it++) {
        if ((*it).first.substr(0,l) != kwex) continue;
        if (fgd && !isdigit((*it).first[l])) continue;
        lstsup.push_back((*it).first);
      }
      list<string>::iterator it2;
      for (it2 = lstsup.begin(); it2 != lstsup.end(); it2++)
        dvl.DeleteKey(*it2);
    }
  }
  dvl_fitskw->Merge(dvl);
}

/* Nouvelle-Fonction */
void SKM_MergeFITSKeywords2(DVList & dvl)
{
  // Cleaning the keywords 
#define SZexlst 21
  const char *exlst[SZexlst]=
  {"SIMPLE","BITPIX" ,"NAXIS" ,"NAXIS#" ,"PCOUNT","GCOUNT",
     "EXTEND","ORIGIN" ,"DATE*" ,"TFIELDS","TTYPE#","TFORM#",
     "TUNIT#","EXTNAME","CTYPE#","CRVAL#" ,"CRPIX#","CDELT#",
     "XTENSION","INSTRUME","TELESCOP"};
  char  kwex[32];
  int i,l;
  for (i=0; i<SZexlst; i++) {
    strncpy(kwex, exlst[i], 32);
    l = strlen(kwex)-1;
    if ((kwex[l] != '*') && (kwex[l] != '#')) {
      dvl.DeleteKey(kwex);
    } 
    else {
      bool fgd = (kwex[l] == '#') ? true : false;
      list<string> lstsup;
      kwex[l] = '\0';
      DVList::ValList::const_iterator it;
      for (it = dvl.Begin(); it != dvl.End(); it++) {
        if ((*it).first.substr(0,l) != kwex) continue;
        if (fgd && !isdigit((*it).first[l])) continue;
        lstsup.push_back((*it).first);
      }
      list<string>::iterator it2;
      for (it2 = lstsup.begin(); it2 != lstsup.end(); it2++)
        dvl.DeleteKey(*it2);
    }
  }
  dvl_fitskw->Merge(dvl);
}

/* Nouvelle-Fonction */
void RadSpec2Nt(RadSpectra & rs, POutPersist & so, string name)
{
  const char *ntn[2] = {"nu","fnu"};
  NTuple nt(2,ntn);  // Creation NTuple (AVEC new )
  float xnt[2];
  double nu;
  double numin = rs.minFreq();
  double numax = rs.maxFreq();
  int nmax = 500;
  double dnu = (numax-numin)/nmax;
  for(int k=0; k<nmax; k++) {
    nu = numin+k*dnu;
    xnt[0] = nu;
    xnt[1] = rs.flux(nu);
    nt.Fill(xnt);
  }
  ObjFileIO<NTuple> oiont(&nt);
  oiont.Write(so, name);
  return;
}

/* Nouvelle-Fonction */
void SpectralResponse2Nt(SpectralResponse& sr, POutPersist & so, string name)
{
  const char *ntn[2] = {"nu","tnu"};
  NTuple nt(2,ntn);  // Creation NTuple (AVEC new )
  float xnt[2];
  double nu;
  double numin = sr.minFreq();
  double numax = sr.maxFreq();
  int nmax = 500;
  double dnu = (numax-numin)/nmax;
  for(int k=0; k<nmax; k++) {
    nu = numin+k*dnu;
    xnt[0] = nu;
    xnt[1] = sr.transmission(nu);
    nt.Fill(xnt);
  }
  ObjFileIO<NTuple> oiont(&nt);
  oiont.Write(so, name);
  return;
}