#include "machdefs.h"
#include <stdlib.h>
#include <stdio.h>

#include <iostream.h>
#include <math.h>

#include <string>
#include <vector>

#include "timing.h"
#include "sambainit.h"
#include "pexceptions.h"
#include "datacards.h"
#include "fitsioserver.h"

#include "spherehealpix.h"
#include "fiospherehealpix.h"

#include "radspecvector.h"
#include "blackbody.h"
#include "derivblackbody.h"
#include "nupower.h"

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

// -----------------------------------------------------------------
// ------------- Function declaration ------------------------------
int CheckCards(DataCards & dc, string & msg);
char * BuildFITSFileName(string const & fname);
SpectralResponse * getSpectralResponse(DataCards & dc);
RadSpectra * getEmissionSpectra(DataCards & dc, int nc);
void RadSpec2Nt(RadSpectra & rs, POutPersist & so, string name);
void SpectralResponse2Nt(SpectralResponse& sr, POutPersist & so, string name);

// 
// treating maps
//---------------------------------------------------------
template <class T> void  MeanSig(NDataBlock<T> const & dbl, double& gmoy, double& gsig);
template <class T> void Sph2Sph(PixelMap<T>& in, PixelMap<T>& out);
//
// 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 integratedMap(SpectralResponse&  sr,   
  //			      PixelMap<T>& betaMap, 
  //			      double normFreq, 
  //			      PixelMap<T>& outMap);
template <class T>
void addComponentBeta(SphereHEALPix<T>& finalMap, 
		      SphereHEALPix<T>& mapToAdd,
		      SphereHEALPix<T>& intBetaMap, double K);
//
// -----------------------------------------------------------------

//  ----- 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: skymixer parameterFile outputfitsname [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 << " !!!! 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");

SphereHEALPix<float> outgs(hp_nside);
bool okout = false;

try {
  if (rdmap) {  // Reading map from FITS file
    FitsIoServer fios;
    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;
    fios.load(outgs, ifnm, 2);
    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 {
    cout << " Output HealPix Map  created - NbPixels= " << outgs.NbPixels() << endl;
    outgs.SetPixels(0.);
  }

  // Decoding detection pass-band filter 
  sr = getSpectralResponse(dc);
  PrtTim(" After FilterCreation ");

  FitsIoServer fios; // Our FITS IO Server
  char * flnm, buff[90];
  string key;

  SphereHEALPix<float> ings(hp_nside);  // The input map
  double K = 1.;
  double freqOfMap = 1.;
  // Loop over sky component 
  int sk;
  for(sk = 0; sk<nskycomp; sk++) {
    cout << " Processing sky component No " << sk+1 << endl;
    
    sprintf(buff, "%d", sk+1);
    key = (string)"MAPFITSFILE" + buff;
    flnm = BuildFITSFileName(dc.SParam(key, 0));
    
    K = dc.DParam(key, 1, 1.);
    
    bool mapDependentOfFreq = false;
    key = (string)"BETAFITSFILE"+ buff;
    if(dc.HasKey(key)) 
      {
	mapDependentOfFreq = true;
      }
    cout << "the map has a mapDependentOfFreq = " << 
      mapDependentOfFreq << endl;

    cout << " Reading Input FITS map " << (string)flnm << endl;
    fios.load(ings, flnm, 2);
    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;
    }
    // getting Emission spectra  
    if(!mapDependentOfFreq)
      {
	if (es) { delete es; es = NULL; }
	es = getEmissionSpectra(dc, sk);
	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;
	fios.load(betaMap, flnm, 2);
	if(nSideForInt<0) nSideForInt = sqrt(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);
	    betaMap.Resize(8);
	    MeanSig(intBetaMap.DataBlock(), moy, sig );
	    if (printlev > 4) 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;
      }

    okout = true;
    if (printlev > 1) {
      MeanSig(outgs.DataBlock(), moy, sig );
      cout << " MeanSig for Sum map - Mean= " << moy << " Sigma= " << sig << endl;
    }
    sprintf(buff, "End of Proc. Comp. %d ", sk+1);
    cout << " --------------------------------------- \n" << endl;    
    PrtTim(buff);
  }   // End of sky component loop 
}   // End of try block

catch (PException exc) {
  msg = exc.Msg();
  cerr << " !!!! skymixer - Catched exception - Msg= " << exc.Msg() << endl;
  rc = 50;
  }   

// Saving the output map in FITS format 
if (okout) {
  FitsIoServer fios;
  fios.save(outgs, arg[2]);
  cout << "Output Map (SphereHEALPix<float>) written to FITS file " 
       << (string)(arg[2]) << endl;
  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 ");
   }
  }
  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';
hp_nside = 32;
nskycomp = 0;
debuglev  = 0;
printlev = 0;

int rc = 0;
string key, key2;

  //  Cheking datacards
  if (dc.NbParam("SKYMIX") < 2)   {
     rc = 71; 
     msg = "Invalid parameters  - Check @SKYMIX card ";
     return(rc);
  }
  key = "READMAP";
  if (dc.HasKey(key)) {
    if (dc.NbParam(key) < 1) {
     rc = 72; 
     msg = "Invalid parameters  - Check @READMAP card ";
     return(rc);
    }
   else rdmap = true;
  }

//  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); 
    }

  // 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 key3;
  string key4;
  string key5;
  for(kc=0; kc<ncomp; kc++) {
    sprintf(buff, "MAPFITSFILE%d", kc+1);
    key = buff;
    if (dc.NbParam(key) < 1 ) {   
      msg = "Missing or invalid card : " + key + " " + key2 ;
      pb = true;  break;
    }
    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;
    if ( (dc.NbParam(key) < 3) && (dc.NbParam(key2) < 1) && (dc.NbParam(key3) < 6) && (dc.NbParam(key4)<2)
	 && (dc.NbParam(key5)<1))  {   
      msg = "Missing card or invalid parameters : " + key + " " + key2 + " " + key3 + " " + key4+ " " + key5;
      pb = true;  break;
    }

    }

  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)
{
  SpectralResponse * filt = NULL;

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

  if (dc.HasKey(key) ) {      // Reading FITS filter file
    FitsIoServer fios;
    char ifnm[256];
    strncpy(ifnm, dc.SParam(key, 1).c_str(), 255);
    ifnm[255] = '\0';
    Matrix mtx(2,10);
    fios.load(mtx, ifnm);
    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);
}

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

  string key = (string)"SPECTRAFITSFILE" + numb;
  string key2 = (string)"BLACKBODY" + numb;
  string key5 = (string)"DIPOLE" + 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;
    FitsIoServer fios;
    Matrix mtx(2,10);
    fios.load(mtx, ifnm);
    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, 2.726));
    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);
}




/* 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(sqrt(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)
{
  //  int nbpix = intBetaMap.NbPixels();
  //  SphereHEALPix<T> newMap(sqrt(nbpix/12));
  Sph2Sph(betaMap,intBetaMap);
  for(int i=0; i<intBetaMap.NbPixels(); i++)
    {
      RadSpectra* rs =  new PowerLawSpectra
	(1.,-intBetaMap(i), 0., normFreq);
      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())
    {
      cout << "nbpix != intBetaMap.NbPixels()" << endl;
      SphereHEALPix<T> bigBetaMap(sqrt(nbpix/12));
      cout << "new map with size corresponding to mapToAdd" <<  sqrt(nbpix/12) << endl;
      Sph2Sph(intBetaMap,bigBetaMap);
      for(int i=0; i<finalMap.NbPixels();i++)
 	{
 	  finalMap(i) += coeff*mapToAdd(i)*bigBetaMap(i);
 	}
    }
  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 */

template <class T>
void Sph2Sph(PixelMap<T>& in, PixelMap<T>& out)
  // Cette fonction remplit la sphere out a partir de la sphere in
  // Les spheres peuvent etre de type et de pixelisations differentes
{
  int k,kin,kout;
  double teta,phi;
  
  int* cnt = new int[out.NbPixels()+1];
  for(kout=0; kout<out.NbPixels(); kout++)
    { cnt[kout] = 0; out(kout) = 0.; }
  
  for(kin=0; kin<in.NbPixels(); kin++) {
    in.PixThetaPhi(kin, teta, phi);
    kout = out.PixIndexSph(teta, phi);
    out(kout) += in(kin);
    
    cnt[kout] ++;
    
  }
  
  double moy, sig, dcn;
  moy = 0.; sig = 0.;
  for(kout=0; kout<out.NbPixels(); kout++) {
    dcn = cnt[kout];  moy += dcn;	 sig += (dcn*dcn);
    if (cnt[kout] > 0)  out(kout) /= dcn;
    else {
      out.PixThetaPhi(kout, teta, phi);
      int pixel = in.PixIndexSph(teta,phi);
      out(kout) = in.PixVal(pixel);
    }
  }
  
  moy /= out.NbPixels();
  sig = sig/out.NbPixels() - moy*moy;
  if (sig >= 0.)	sig = sqrt(sig);
  
  delete[] cnt;
  
  printf("Sph2Sph_Info NbPix In= %d  Out= %d CntMoy,Sig= %g %g\n",
	 in.NbPixels(), out.NbPixels(), moy, sig);
  PrtTim("End of  Sph2Sph() ");
}



/* 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 */
void RadSpec2Nt(RadSpectra & rs, POutPersist & so, string name)
{
  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)
{
  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;
}