Merge branch 'u/jct/dev' into '48-improve-output-key-management-and-header-construction'

# Conflicts:
#   source/photoz_lib.cpp
#   tests/test_suite.sh

source/PDF.cpp

@@ -7,14 +7,13 @@
 
 #include "globals.h"
 #include "PDF.h"
+#include "cosmology.h"
 #include <fstream> // print output file
 #include <iostream> // print standard file
 #include <sstream>
 #include <string>
 #include <vector>
 #include <cmath>
-#include <algorithm>    // sort
-
 
 using namespace std;
 
@@ -23,7 +22,24 @@ using namespace std;
 /*
  CONSTRUCTOR OF THE PDF
 */
-void PDF::initialize(vector<double> xaxisIn) {
+PDF::PDF(const double min, const double step, const size_t size): PDF() {
+  //by construction
+  scaleLinear = 1; 
+  scaleStep = step;
+  scaleMin = min;
+  vPDF.resize(size, 0.0);
+  chi2.resize(size, 1.e9);
+  ind.resize(size, 0);
+  xaxis.resize(size, 0.);
+  for(size_t k=0; k<size; k++) {
+    xaxis[k] = min + double(k)*step;
+  }
+  scaleMax = xaxis.back();
+  vsize = xaxis.size();
+  return;
+}
+
+void PDF::initialize(const vector<double>& xaxisIn) {
 
   // define the size of the PDF vector
   vsize=int(xaxisIn.size());
@@ -71,7 +87,6 @@ void PDF::normalization() {
     for (size_t k=0; k<vPDF.size()-1; k++){
       somme+=(vPDF[k]+vPDF[k+1])/2.*(xaxis[k+1]-xaxis[k]);
     }
-
     // Renormalize the PDF
     if(somme>0){
       for(int k=0;k<vsize;k++){vPDF[k]=vPDF[k]/somme; }
@@ -98,7 +113,7 @@ void PDF::chi2toPDF() {
 /*
  INDEX OF THE MINIMUM CHI2 VALUE
 */
-int PDF::chi2mini() {
+int PDF::chi2mini() const {
 
 
   double chi2min=1.e50;
@@ -123,7 +138,7 @@ int PDF::chi2mini() {
  Description : Parabolic interpolation for minmum Chi2
  Origine     : Bevington book
 */
-double PDF::int_parab() {
+double PDF::int_parab() const {
 
   double xmin=-99;
   double dxb,dxa,num,den;
@@ -167,7 +182,7 @@ double PDF::int_parab() {
 /*
  Find the values of xaxis crossing chiMin+dchi
 */
-vector<double> PDF::uncMin(double dchi){
+vector<double> PDF::uncMin(double dchi) const {
 
   double chiLim; // chi2 limit to be considered
 
@@ -211,12 +226,61 @@ vector<double> PDF::uncMin(double dchi){
 }
 
 
-
+pair<double, double> PDF::credible_interval(float level){
+  pair<double, double> result;
+  double bound_left_id, bound_right_id;
+  double x1, x2, y1, y2,y;
+    //version without any interpolation for now
+  if(level>1.) level /= 100.;
+  //the logic of the algorithm is the following
+  //Find Mode
+  // - Find Mode
+  size_t maxid = get_maxid();
+  double xmax = xaxis[maxid];
+  // - Compute the full cumulant
+  vector<double> cumulant;
+  double tmp=0;
+  cumulant.push_back(0);
+  for(size_t k=0; k<xaxis.size()-1;k++){
+    tmp += (xaxis[k+1]-xaxis[k])*(vPDF[k+1]+vPDF[k])/2.;
+    cumulant.push_back(tmp);
+  }
+  if((cumulant.back()-1.)> 1.e-6){
+    for(auto& c : cumulant) c /= cumulant.back();
+  }
+  double cumul_left = cumulant[maxid];
+  // - Use cumulant to compute integral from Mode to end
+  double cumul_right = 1.-cumul_left;
+  // - If higher than xval/2, then interpolate cumulant to find
+  // where delta cumulant is level/2, and set new_val to level/2
+  if(cumul_right>level/2.){
+    bound_right_id =
+      upper_bound(cumulant.begin(), cumulant.end(), cumul_left+level/2.) - cumulant.begin();
+    // args in reverse order y, y1, x1, y2, x2, in order to get x instead of y
+    result.second = linear_interp(cumul_left+level/2., cumulant[maxid],
+				  xmax, cumulant[bound_right_id], xaxis[bound_right_id]);
+    bound_left_id =
+      lower_bound(cumulant.begin(), cumulant.end(), cumul_left-level/2.) - cumulant.begin()-1;
+    result.first = linear_interp(cumul_left-level/2., cumulant[bound_left_id] ,xaxis[bound_left_id],
+				 cumulant[maxid], xmax);
+  }
+  // - Else compute delta integrant from Mode to end, set as high_bound return value,
+  // and set new_val to level - high_bound
+  else {
+    result.second = xaxis.back();
+    bound_left_id =
+      lower_bound(cumulant.begin(), cumulant.end(), 1.-level) - cumulant.begin()-1;
+    result.first = xaxis[bound_left_id];
+    result.first = linear_interp(1-level, cumulant[bound_left_id] ,xaxis[bound_left_id],
+				 cumulant[maxid], xmax); 
+  }
+  return result;
+}
   
 /*
  Find the values of xaxis including xvalue% of the PDF area
 */
-double PDF::probaBay(double xvalue){
+double PDF::probaBay(double xvalue) const {
 
 
   double result;
@@ -268,43 +332,6 @@ double PDF::probaBay(double xvalue){
 
 }
 
-//
-// Find the index in the PDF of the considered value
-// If out of range, output a negative value 
-int PDF::index(const double inVal){
-
-  int res=-1;
-
-  // If the value is well with the allowed scale
-  if(inVal>=scaleMin && inVal<=scaleMax){
-
-   // If the scale is linear, the index is computed easily
-   if(scaleLinear == 1){
-       res=int(round((inVal-scaleMin)/scaleStep));
-   }else{
-
-    // If the scale is not linear, go step by step
-    // Loop on the xaxis grid
-    for(int k=0; k<vsize;k++){
-     // if the input value match the one of the grid
-     if (inVal>=xaxis[k] && inVal<xaxis[k+1]){
-      // result is found
-      res=k;
-      break;
-     }
-    }       
-   }
-
-  }else{
-    //cout << " Value is out of range for the PDF:" << inVal << " " << scaleMin  << " " <<scaleMax<< endl;
-    if(inVal<scaleMin){res=0;}else{res=vsize-1;}
-  }
-
-  return res;
-
-}
-  
-
 
 //
 // Search for high peaks in The ML function vs xaxis

source/PDF.h

@@ -9,6 +9,9 @@
 
 #include <string>
 #include <vector>
+#include <stdexcept>
+#include <algorithm>    // sort
+#include "cosmology.h"
 
 using namespace std;
 
@@ -42,18 +45,47 @@ class PDF
 
 
   //Prototype
+  PDF(const double min, const double step, const size_t size); ///< initialize the pdf with minimal value, step, and total number of bins.
   void normalization(); ///< normalize the pdf to 1.
   void chi2toPDF(); ///< compute \f$\exp(-\chi^2/2)\f$ from the vector of \f$\chi^2\f$ values
-  int chi2mini(); ///< index of the minimum \f$\chi^2\f$ value
-  double int_parab(); ///< compute the local \f$\chi^2\f$ by parabolic interpolation
-  vector<double> uncMin(double dchi); ///< compute the bounds of the confidence interval defined by \f$\chi^2_{min} + dchi\f$
-  double probaBay(double xvalue); ///< return the inverse cumulative function at xvalue
-  int index(const double inVal); ///< return the vector index where the pdf vector is closest to inVal
-  void initialize(vector<double> xaxisIn); ///< initialize the pdf representation over the redshift grid xaxisIn
+  int chi2mini() const; ///< index of the minimum \f$\chi^2\f$ value
+  double int_parab() const; ///< compute the local \f$\chi^2\f$ by parabolic interpolation
+  vector<double> uncMin(double dchi) const; ///< compute the bounds of the confidence interval defined by \f$\chi^2_{min} + dchi\f$
+  double probaBay(double xvalue) const; ///< return the inverse cumulative function at xvalue
+  void initialize(const vector<double>& xaxisIn); ///< initialize the pdf representation over the redshift grid xaxisIn
   void secondMax(const double win); ///< search for high peaks in The ML function vs xaxis, and sort them from highest to smallest peaks in ML
 
+  /// return the vector index where the pdf vector is closest to inVal
+  inline int index(const double inVal){
+    if(size()==0)
+      throw invalid_argument("PDF::index error : xaxis has size 0" );
+    return indexz(inVal, xaxis);
+  } 
+  inline void setVal(const vector<double>& vals) {for(size_t i=0;i<vals.size();i++) vPDF[i]=vals[i];}
 
-};
+  inline size_t size() const {return vsize;}
+  
+  inline double getX(const size_t index) const {
+    if(index>=size()) throw invalid_argument("PDF::getX : index >= size()" );
+    return xaxis[index];
+  }
+  
+  inline double getY(const size_t index) const {
+    if(index>=size()) throw invalid_argument("PDF::getY : index >= size()" );
+    return vPDF[index];
+  }
 
+  inline size_t get_maxid() const {return max_element(vPDF.begin(), vPDF.end())-vPDF.begin();}
+  inline double get_max() const {return xaxis[get_maxid()];}
+  inline double linear_interp(double x, double x1, double y1, double x2, double y2){
+    return y1 + (x-x1)*(y2-y1)/(x2-x1);
+  }
+  //for marginalized pdf use
+  double int_parab2() {for(size_t k=0; k<vPDF.size(); k++){chi2[k]=-vPDF[k];} return int_parab();}
+
+  pair<double, double> credible_interval(float xval);
+
+
+};
 
 #endif

source/SED.h

 /*
- 17/11/2014
- Class to store one SED
+  17/11/2014
+  Class to store one SED
 */
 
 // avoid multiple def of the same class
@@ -23,105 +23,105 @@
  */
 class SED{
 protected:
-    int idAge; ///< index of the age for this SED object
-    string type; ///< type of the SED object: S(tar), G(alaxy), Q(SO) (lower case allowed)
+  int idAge; ///< index of the age for this SED object
+  string type; ///< type of the SED object: S(tar), G(alaxy), Q(SO) (lower case allowed)
 
 public:
-    vector<oneElLambda> lamb_flux; ///< vector of oneElLambda elements, capturing the SED template
-    vector<double> kcorr, ///< k-correction term
-	kcolor, ///< k-correction + color rest-frame term (to minimize the SED dependency)
-	mag_z0; ///< magnitude at z=0
-	vector<double> mag,flux, ///< magnitude and flux of the model
-	flux_ori, ///< original flux vector before any transformation is applied
-	mag_ori; ///< original mag vector before any transformation is applied
-	string name;
+  vector<oneElLambda> lamb_flux; ///< vector of oneElLambda elements, capturing the SED template
+  vector<double> kcorr, ///< k-correction term
+    kcolor, ///< k-correction + color rest-frame term (to minimize the SED dependency)
+    mag_z0; ///< magnitude at z=0
+  vector<double> mag,flux, ///< magnitude and flux of the model
+    flux_ori, ///< original flux vector before any transformation is applied
+    mag_ori; ///< original mag vector before any transformation is applied
+  string name;
   bool has_emlines; ///< True if the emission lines have been computed, false if not
-	int nummod;
-        int nlib,index,index_z0,nbusIR;
+  int nummod;
+  int nlib,index,index_z0,nbusIR;
   double red,chi2,prob,dm,lnir,luv,lopt;
-	double mass,age,sfr,ltir; // need to put it out of GalSED since used in the PDF without knowing that it's a gal.
-	double ebv,mag0;
-	int extlawId;
-	double qi[4]; ///< Store the number flux (phot/cm\f$^{-2}\f$s\f$^{-1}\f$) of ionizing photons for HeII, HeI, H, and H2. See SED::calc_ph. In practice, qi[2] only is used, and only for the physical modeling of emission lines (EM_LINES="PHYS", see GalMag::read_SED)
-        vector <oneElLambda> fac_line,fac_line_ori;
+  double mass,age,sfr,ltir; // need to put it out of GalSED since used in the PDF without knowing that it's a gal.
+  double ebv,mag0;
+  int extlawId;
+  double qi[4]; ///< Store the number flux (phot/cm\f$^{-2}\f$s\f$^{-1}\f$) of ionizing photons for HeII, HeI, H, and H2. See SED::calc_ph. In practice, qi[2] only is used, and only for the physical modeling of emission lines (EM_LINES="PHYS", see GalMag::read_SED)
+  vector <oneElLambda> fac_line,fac_line_ori;
   
-        // Constructors defined in SED.cpp
-	SED(const string nameC, int nummodC=0, string typeC="N");
-	SED(const string nameC, double tauC, double ageC, int keepAgeC, string formatC, int nummodC, string typeC, int idAgeC);
-	SED(SED const& p){
-		idAge = p.idAge;
-		type = p.type;
-		lamb_flux = p.lamb_flux;
-		kcorr = p.kcorr;
-		kcolor = p.kcolor;
-		mag = p.mag;
-		flux = p.flux;
-		flux_ori = p.flux_ori;
-		mag_ori = p.mag_ori;
-		name = p.name;
-		has_emlines = p.has_emlines;
-		nummod = p.nummod;
-		nlib = p.nlib;
-		index = p.index;
-		index_z0 = p.index_z0;
-		red = p.red;
-		chi2 = p.chi2;
-		dm = p.dm;
-		mass = p.mass;
-		age = p.age;
-		sfr = p.sfr;
-		ltir = p.ltir;
-		ebv = p.ebv;
-		extlawId = p.extlawId;
-                fac_line=p.fac_line;
-	};
-	virtual ~SED();
-    ///\brief Read fltFile assumed to be ASCII and build the #lamb_flux vector of oneElLambda elements.
-    ///
-    /// Negative flux values are set to 0, oneElLambda::ori=1 to indicate that it is a SED.
-    /// The #lamb_flux vector is filled iteratively with each line of the file, and is finally
-    /// sorted by ascending lambda.
-    /// More complex input types are treated in inherited class methods.
-    void read(string fltFile);
-    void warning_integrateSED(const flt& filter);
-    vector<double>  integrateSED(const flt& filter);
-    void resample(vector<oneElLambda>& lamb_all, vector<oneElLambda>& lamb_new, const int origine, const double lmin, const double lmax);
+  // Constructors defined in SED.cpp
+  SED(const string nameC, int nummodC=0, string typeC="N");
+  SED(const string nameC, double tauC, double ageC, int keepAgeC, string formatC, int nummodC, string typeC, int idAgeC);
+  SED(SED const& p){
+    idAge = p.idAge;
+    type = p.type;
+    lamb_flux = p.lamb_flux;
+    kcorr = p.kcorr;
+    kcolor = p.kcolor;
+    mag = p.mag;
+    flux = p.flux;
+    flux_ori = p.flux_ori;
+    mag_ori = p.mag_ori;
+    name = p.name;
+    has_emlines = p.has_emlines;
+    nummod = p.nummod;
+    nlib = p.nlib;
+    index = p.index;
+    index_z0 = p.index_z0;
+    red = p.red;
+    chi2 = p.chi2;
+    dm = p.dm;
+    mass = p.mass;
+    age = p.age;
+    sfr = p.sfr;
+    ltir = p.ltir;
+    ebv = p.ebv;
+    extlawId = p.extlawId;
+    fac_line=p.fac_line;
+  };
+  virtual ~SED();
+  ///\brief Read fltFile assumed to be ASCII and build the #lamb_flux vector of oneElLambda elements.
+  ///
+  /// Negative flux values are set to 0, oneElLambda::ori=1 to indicate that it is a SED.
+  /// The #lamb_flux vector is filled iteratively with each line of the file, and is finally
+  /// sorted by ascending lambda.
+  /// More complex input types are treated in inherited class methods.
+  void read(string fltFile);
+  void warning_integrateSED(const flt& filter);
+  vector<double>  integrateSED(const flt& filter);
+  void resample(vector<oneElLambda>& lamb_all, vector<oneElLambda>& lamb_new, const int origine, const double lmin, const double lmax);
   
-    ///\brief Generate a calibration SED based on the argument calib
-    ///
-    ///@param lmin start of the lambda vector
-    ///@param lmax end of the lambda vector
-    ///@param Nsteps number of intervals between $lambda values (hence there are Nsteps+1 values of \lambda)
-    ///@param calib: parameter FILTER_CALIB passed as argument to define the calibration function \f$C(\lambda)\f$
-    /// - calib=0 : \f$C(\lambda)=\lambda^{-2}\f$
-    /// - calib=1 : \f$C(\lambda)=\lambda^{-1}\f$
-    /// - calib=2 : \f$C(\lambda)=\lambda^{-3}\f$
-    /// - calib=3 : \f$C(\lambda)=Blackbody(\lambda, T=10000K)\f$
-    /// - calib=4 : \f$C(\lambda)=Blackbody(\lambda, T=10000K)\f$
-    /// - calib=5 : \f$C(\lambda)=\lambda^{-3}\f$
-    void generateCalib(double lmin, double lmax, int Nsteps, int calib);
-    /// return the size of the internal vector #lamb_flux
-    int nlines();
-    /// rescale the lamb_flux.val as val *= scaleFac
-    void rescale(double scaleFac);
-    ///apply  \a #shifts to the magnitude : mag += shifts
-    void applyShift(const vector<double>& shifts, const int imagm);
-    /// compute and return \f$f = 10^{-0.4(mag+48.6)}\f$  
-    void mag2flux();
-    /// Store the original flux and mag into  \a #flux_ori and \a #mag_ori, before any transformation (like emission line addition) is applied. 
-    void keepOri();
-    void writeType(){cout << "SED Type " << type << endl;};
-
-    void sumSpectra(SED addSED);
+  ///\brief Generate a calibration SED based on the argument calib
+  ///
+  ///@param lmin start of the lambda vector
+  ///@param lmax end of the lambda vector
+  ///@param Nsteps number of intervals between $lambda values (hence there are Nsteps+1 values of \lambda)
+  ///@param calib: parameter FILTER_CALIB passed as argument to define the calibration function \f$C(\lambda)\f$
+  /// - calib=0 : \f$C(\lambda)=\lambda^{-2}\f$
+  /// - calib=1 : \f$C(\lambda)=\lambda^{-1}\f$
+  /// - calib=2 : \f$C(\lambda)=\lambda^{-3}\f$
+  /// - calib=3 : \f$C(\lambda)=Blackbody(\lambda, T=10000K)\f$
+  /// - calib=4 : \f$C(\lambda)=Blackbody(\lambda, T=10000K)\f$
+  /// - calib=5 : \f$C(\lambda)=\lambda^{-3}\f$
+  void generateCalib(double lmin, double lmax, int Nsteps, int calib);
+  /// return the size of the internal vector #lamb_flux
+  int nlines();
+  /// rescale the lamb_flux.val as val *= scaleFac
+  void rescale(double scaleFac);
+  ///apply  \a #shifts to the magnitude : mag += shifts
+  void applyShift(const vector<double>& shifts, const int imagm);
+  /// compute and return \f$f = 10^{-0.4(mag+48.6)}\f$  
+  void mag2flux();
+  /// Store the original flux and mag into  \a #flux_ori and \a #mag_ori, before any transformation (like emission line addition) is applied. 
+  void keepOri();
+  void writeType(){cout << "SED Type " << type << endl;};
+
+  void sumSpectra(SED addSED);
 
   void reduce_memory(vector<flt> allFlt);
   
   /*
-     * These functions are different depending on the type of SED
-     */
-    virtual void fit(const onesource & source,const int imagm, const double lmasi, const double lmass,const double magabsB[2], const double magabsF[2], const double funz0, const int bp[2]){};
-    virtual void fitIR(const onesource& source,const int imagm, const double consideredZ, const string fit_frsc, cosmo lcdm){};
-    virtual void writeSED(ofstream& ofs, ofstream& ofsPhys, ofstream& ofsDoc);
+   * These functions are different depending on the type of SED
+   */
+  virtual void fit(const onesource & source,const int imagm, const double lmasi, const double lmass,const double magabsB[2], const double magabsF[2], const double funz0, const int bp[2]){};
+  virtual void fitIR(const onesource& source,const int imagm, const double consideredZ, const bool fit_frsc, cosmo lcdm){};
+  virtual void writeSED(ofstream& ofs, ofstream& ofsPhys, ofstream& ofsDoc);
   inline void writeSED(const string& binFile, const string& physFile, const string& docFile) {
     ofstream sdocOut, sphysOut, sbinOut;
     sdocOut.open(docFile.c_str());
@@ -148,14 +148,14 @@ public:
     if(!sbinIn){throw invalid_argument("Can't open file " + fname); }
     readSEDBin(sbinIn);
   }
-    /// read the SED library when it is in binary format
-    virtual void readSEDBin(ifstream& ins);
-    virtual void readMagBin(ifstream& ins){};
-    virtual void sumEmLines(){};
-    /// for each magnitude \a #mag[k] compute kcorr = mag[k] - mag_z0[k] - distMod
-    virtual void kcorrec(){};
-    virtual void prop(){};
-    virtual void add_neb_cont(){}; // Add continuum
+  /// read the SED library when it is in binary format
+  virtual void readSEDBin(ifstream& ins);
+  virtual void readMagBin(ifstream& ins){};
+  virtual void sumEmLines(){};
+  /// for each magnitude \a #mag[k] compute kcorr = mag[k] - mag_z0[k] - distMod
+  virtual void kcorrec(){};
+  virtual void prop(){};
+  virtual void add_neb_cont(){}; // Add continuum
   /*!
    * Compute the number flux of photons able to ionize HeII, HeI, H, and H2
    * For a given SED, this amounts to compute the integral
@@ -167,17 +167,17 @@ public:
    *
    * Results are stored in the q_i array member of size 4 of the SED instance.
    */
-    virtual void calc_ph(){}; // Number of inoizing photons
-    virtual void SEDproperties(){};
-
-    //in read_lib for zphota
-    virtual void clean(){lamb_flux.clear();mag.clear();mag_ori.clear();flux_ori.clear();flux.clear();kcorr.clear();kcolor.clear();fac_line.clear();};
-    virtual void setOthers(){};
-    void fit_normalization(const onesource& source, const int imagm);
-    void fit_normalizationIR(const onesource& source, const int imagm, const double dmcor);
-    void fit_chi2(const onesource& source, const int imagm);
-    void fit_chi2IR(const onesource& source, const int imagm, const string fit_frsc, const double dmcor);
-    inline bool is_same_model(const SED& other) {return ((*this).nummod == other.nummod && (*this).ebv == other.ebv && (*this).age == other.age) ;}
+  virtual void calc_ph(){}; // Number of inoizing photons
+  virtual void SEDproperties(){};
+
+  //in read_lib for zphota
+  virtual void clean(){lamb_flux.clear();mag.clear();mag_ori.clear();flux_ori.clear();flux.clear();kcorr.clear();kcolor.clear();fac_line.clear();};
+  virtual void setOthers(){};
+  void fit_normalization(const onesource& source, const int imagm);
+  void fit_normalizationIR(const onesource& source, const int imagm, const double dmcor);
+  void fit_chi2(const onesource& source, const int imagm);
+  void fit_chi2IR(const onesource& source, const int imagm, const bool fit_frsc, const double dmcor);
+  inline bool is_same_model(const SED& other) {return ((*this).nummod == other.nummod && (*this).ebv == other.ebv && (*this).age == other.age) ;}
 };
 
 
@@ -187,77 +187,77 @@ public:
  */
 class ExtSED : public SED{
 public:
-	double distMod;
-        int extNum;
+  double distMod;
+  int extNum;
     
-       // Constructors defined in SED.cpp
-	ExtSED(const string nameC, int nummodC=0, string typeC="N");
-	ExtSED(const string nameC, double tauC, double ageC, int keepAgeC, string formatC, int nummodC, string typeC, int idAgeC);
+  // Constructors defined in SED.cpp
+  ExtSED(const string nameC, int nummodC=0, string typeC="N");
+  ExtSED(const string nameC, double tauC, double ageC, int keepAgeC, string formatC, int nummodC, string typeC, int idAgeC);
 
-        ExtSED(SED const& p):SED(p) {}; // this is dangerous to have around!!
-        ExtSED(ExtSED const& p):SED(p),distMod(p.distMod) {};
-	~ExtSED();
+  ExtSED(SED const& p):SED(p) {}; // this is dangerous to have around!!
+  ExtSED(ExtSED const& p):SED(p),distMod(p.distMod) {};
+  ~ExtSED();
 
-	void redshift();
-	void applyExt(const ext& oneext);
-	void applyExtLines(const ext& oneext);
-	void applyOpa(const vector<opa>& opaAll);
+  void redshift();
+  void applyExt(const ext& oneext);
+  void applyExtLines(const ext& oneext);
+  void applyOpa(const vector<opa>& opaAll);
 
-	void clean(){SED::clean(); };
+  void clean(){SED::clean(); };
 };
 
 /// concrete ExtSED implementation for galaxy objects
 class GalSED : public ExtSED{
 public:
-        vector<double> flEm;
- 	string format;
-	double tau,zmet,d4000,fracEm;
-	int keepAge;
-
-        GalSED(SED const& p):ExtSED(p) {};
-        GalSED(GalSED const& p):ExtSED(p){
-      	    flEm = p.flEm;
-    	    format = p.format;
-    	    tau = p.tau;
-    	    zmet = p.zmet;
-    	    lnir = p.lnir;
-    	    luv = p.luv;
-    	    lopt = p.lopt;
-    	    d4000 = p.d4000;
-    	    fracEm = p.fracEm;
-    	    keepAge = p.keepAge;
-        };
-
-       // Constructors defined in SED.cpp
-	GalSED(const string nameC, int nummodC=0);
-	GalSED(const string nameC, double tauC, double ageC, int keepAgeC, string formatC, int nummodC, string typeC, int idAgeC);
-	~GalSED();
-
-	void SEDproperties();
-	void add_neb_cont();
-	void generateEmEmpUV(double MNUV_int, double NUVR);
-	void generateEmEmpSFR(double MNUV_int, double NUVR);
-	void generateEmPhys(double zmet, double qi); 
-        void generateEmSpectra();
-        void generateEmSpectra(int nstep);
-	void sumEmLines();
-	void kcorrec();
-        void rescaleEmLines();
-        void zdepEmLines(int flag);
-	void prop();
-	void calc_ph();
-
-	void writeSED(ofstream& ofs, ofstream& ofsPhys, ofstream& ofsDoc);
-	void readSEDBin(ifstream& ins);
-
-	void writeMag(bool outasc, ofstream& ofsBin, ofstream& ofsDat, vector<flt> allFilters, string magtyp);
-	void readMagBin(ifstream& ins);
-        void generate_spectra(double zin, double dmin, vector<opa> opaAll);
-
-	void fit(const onesource& source,const int imagm, const double lmasi, const double lmass,const double magabsB[2], const double magabsF[2], const double funz0, const int bp[2]);
-	void fitIR(const onesource& source,const int imagm,  const double consideredZ, const string fit_frsc, cosmo lcdm);
-        double nzprior(const onesource& source,  const int bp[2]);
-	void clean(){ExtSED::clean(); flEm.clear();};
+  vector<double> flEm;
+  string format;
+  double tau,zmet,d4000,fracEm;
+  int keepAge;
+
+  GalSED(SED const& p):ExtSED(p) {};
+  GalSED(GalSED const& p):ExtSED(p){
+    flEm = p.flEm;
+    format = p.format;
+    tau = p.tau;
+    zmet = p.zmet;
+    lnir = p.lnir;
+    luv = p.luv;
+    lopt = p.lopt;
+    d4000 = p.d4000;
+    fracEm = p.fracEm;