photfitter.h

// -*- C++ -*-
// 

#ifndef PHOTFITTER_H
#define PHOTFITTER_H

#include <time.h>

#include <math.h>
#include <assert.h>
#include <fstream>

#include <fileutils.h>
#include "vutils.h"

#include "measuredstar.h"
#include "fittedstar.h"
#include "fittedstarmap.h"
#include "associations.h"



template <class Model>
class PhotFitter 
{
public:
  PhotFitter(Associations& assoc, Model& model);
  
  ~PhotFitter() {}
  
  template<class Op>
  void       LoopOnMeasurements(Op& op) const;
  
  bool         InitFit(int niter=100);
  
  bool         InitModel();
  
  bool         FlipFlopStep();
  
  bool         NewtonRaphsonStep();
  
  bool         Minimize(unsigned int maxiter=100);
  
  int          RemoveOutliers(double nsig);
  
  unsigned int GetNFittedPars() const;
  
  double       GetChi2() const;
  
  int          GetNMes() const;
  
  int          GetNDoF() const;
  
  Vect const&  GetPar() const;
  
  Vect const&  GetParError() const;
  
  Mat const&   GetCovMatrix() const;
  
  Mat const&   GetCorrMatrix() const;
  
  void         DumpNTuples(std::string const& resultdir) const;
  
  void         SetVerbose(bool verbose=true) { verbose_ = verbose; }
  
private:
  bool    ComputeCovMat();
  bool    ComputeCorrMat();
  void    dumpParNTuple_(const std::string& ntname) const;
public:
  void    dumpFittedStarNTuple_(const std::string& ntname) const;
  void    dumpMeasuredStarNTuple_(const std::string& ntname, unsigned int nusrvals=0) const;
private:
  void    updateFittedStars_();
  
  Associations&   assoc_;
  Model&          model_;
  CcdImageList&   imstd::list_;
  FittedStarList& fittedstarstd::list_;
  
  unsigned int npar_;
  Vect         pars_;
  Vect         epars_;
  Mat          covmat_;
  Mat          corrmat_;
  
  double       minstep_;
  
  bool         verbose_;
  mutable bool         cleanChi2_;
  mutable double       chi2_;
  mutable int          nmes_;
  mutable int          ndof_;
};


// 
// template methods
// 

double sq(double x) { return x*x; }

template<class Model>
PhotFitter<Model>::PhotFitter(Associations& assoc, Model& model)
  : assoc_(assoc), model_(model),
    imstd::list_(assoc_.ccdImageList),
    fittedstarstd::list_(assoc_.fittedStarList),
    npar_(0), // one ne sait pas a l'avance combien on aura de parametres...
    minstep_(1.E-3),
    verbose_(false), 
    cleanChi2_(false), 
    chi2_(-1), nmes_(-1), ndof_(-1)
{
  //  cout << " MODEL # pars=" << model_.NPar() << endl;
  //  pars_.allocate(npar_);
  //  epars_.allocate(npar_);
  //  cout << " PHOTFITTER: nfittedstars=" << fittedstarstd::list_.size() << endl;
  
  //  int index = 0;
  //  FittedStarIterator I;
  //  for(I=fittedstarstd::list_.begin();I!=fittedstarstd::list_.end();I++)
  //    {
  //      FittedStar* s = *I;
  //      s->SetIndex(index++);
  //    }
}


template<class Model> template<class Op>
void  PhotFitter<Model>::LoopOnMeasurements(Op& op) const
{
  CcdImageIterator imI;
  MeasuredStarIterator stI;
  
  for(imI=imstd::list_.begin();imI!=imstd::list_.end();imI++)
    {
      CcdImage& ccdImage = **imI;
      MeasuredStarList& catalog = ccdImage.CatalogForFit();
      for(stI=catalog.begin();stI!=catalog.end();stI++)
    {
      MeasuredStar& m = **stI;
      FittedStar*   f = m.GetFittedStar();
      if(!f) continue;
      op(m);
    }
    }
  
  op.Finalize();
}


template<class Model>
class TestModel
{
public:
  TestModel(Model const& model)
    : model_(model), npar_(model_.Size())
  { 
    assert(npar_>0);
    pars_.allocate(npar_);
  }
  
  ~TestModel() { }
  
  void operator()(MeasuredStar const& m)
  {
    model_.CanCompute(m, pars_);
  }
  
  void Finalize() {}

  Vect const GetPars() const { return pars_; }
  
private:
  Model  model_;
  unsigned int npar_;
  Vect pars_;
};



// To initialize the fit, we
// fit separately each group 
// of parameters... If the user
// chose them so that the fit is 
// almost linear, it should give us 
// a good start ...

template<class Model>
bool PhotFitter<Model>::InitFit(int niter) // maybe we should call this FlipFlopFit()
{
  unsigned int i;
  
  //  model_.InitPars(pars_);
  int iter;
  int fcount = 0;
  int rcount = 0;
  npar_ = model_.Size(); 
  //  unsigned int nfittedstars = fittedstarstd::list_.size();
  //  unsigned int nccd = 36;
  //  npar_ = nfittedstars + nccd;
  
  TestModel<Model> tmodel(model_);
  for(iter=0;iter<10;iter++)
    {
      LoopOnMeasurements(tmodel);
      //      rcount = 0; fcount=0;
      //      for(i=0;i<npar_;i++)
      //    {
      //      if(tmodel.GetPars()(i) == 1) rcount++;
      //      if(tmodel.GetPars()(i) == 2) fcount++;
      //    }
      //      cout << " rcount=" << rcount
      //       << " fcount=" << fcount
      //       << endl;
    }
  
  model_.InitIndexTables(tmodel.GetPars());
  
  //  for(i=0;i<npar_;i++)
  //    cout << " pars(" << i << ")=" << tmodel.GetPars()(i) << endl;
  
  //  cout << " model_.NPar() = " << model_.NPar()
  //       << " model_.Npar(0) = " << model_.NPar(0)
  //       << " model_.Npar(1) = " << model_.NPar(1)
  //       << " model_.NPar(2) = " << model_.NPar(2)
  //       << endl;
  
  npar_ = model_.NPar();
  pars_.allocate(npar_);
  model_.InitPars(pars_);
  
  double old_chi2, chi2 = GetChi2();
  int igroup, ngroups = model_.NGroups();
  
  for(iter=0;iter<niter;iter++)
    {
      old_chi2 = chi2;
      FlipFlopStep();
      double chi2 = GetChi2();
      
      // check that the Chi2 is not going up
      if( old_chi2 < chi2 )
    {
      // we should do something clever here
      cout << "[PhotFitter<>::InitFit] Warning increasing chi2 in InitFit() !"
           << endl;
      return false;
    }
      
      // break if too many iterations
      if( iter >= niter )
    {
      cout << "[PhotFitter<>::InitFit] Warning max iter=" << niter
           << " reached" << endl;
      return false;
    }
      
      cout << " Chi2=" << chi2 << endl;
    }
  
  return true;
}



template<class Model>
class MatrixFill
{
public:
  
  MatrixFill(Vect& pars, Mat& M, Vect& B, Model& model, int igroup=-1)
    : npar_( model.NPar(igroup) ),
      npartot_( model.NPar() ),
      igroup_(igroup),
      pars_(pars),
      M_(M), 
      B_(B), 
      model_(model)
  {
    //    assert( pars_.Size() == npar_ ); because pars_.Size() = total # of parameters
    assert( M_.SizeX()   == npar_ );
    assert( M_.SizeY()   == npar_ );
    assert( B_.Size()    == npar_ );
    
    validpar_.allocate(npar_);
    validpar_.Zero();
    deriv_.allocate(npartot_);
    
    //    if(igroup_<0)
    //      {
    //  cout << " igroup=" << igroup_
    //       << " npar_ = " << npar_
    //       << endl;
    //      }
  }
  
  ~MatrixFill() {}
  
  
  void   operator()(MeasuredStar const& m) const
  {
    double val;
    bool status = model_.Compute(m, pars_, val, deriv_);
    if(!status) return;
    
    int g_index = model_.GroupIndex(igroup_);
    
    double flux = m.flux;
    double w    = 1. / sq(m.FluxSig());
    
    unsigned int i, j;
    for(i=0;i<npar_;i++)
      {
    if(deriv_(i+g_index) == 0.) continue;
    B_(i) += -(w * (val-flux) * deriv_(i+g_index));
    for(j=i;j<npar_;j++)
      {
        if(deriv_(j+g_index) == 0.) continue;
        M_(i,j) += w * deriv_(i+g_index) * deriv_(j+g_index);
      }
      }
  }
  
  void   Finalize()
  {
    unsigned int i;
    for(i=0;i<npar_;i++)
      if(M_(i,i) == 0.)
    {
      // we fix the parameter,
      // since it seems to be masked...
      cout << " Finalize: M(" << i << "," << i << ")=" << M_(i,i) << endl;
      M_(i,i) = 1.;
    }
      else
        validpar_(i) = 1;
  }
  
  bool   IsValid(unsigned int i)
  {
    assert(i<npar_);
    return validpar_(i);
    //    return true;
  }
  
private:
  unsigned int npar_;
  unsigned int npartot_;
  int          igroup_;
  
  Vect&  pars_;
  Mat&   M_;
  Vect&  B_;
  Model& model_;
  
  mutable Vect   deriv_;
  Vect   validpar_;
};



template<class Model>
class Chi2Acc 
{
public:
  Chi2Acc(Model& model, Vect const& pars, bool cmpsig=false)
    : chi2_(0.), medres_(0.), sigres_(0.),
      model_(model), pars_(pars), nmes_(0), cmpsig_(cmpsig)
  {
    deriv_.allocate(model_.NPar());
  }
  
  
  ~Chi2Acc() {}
  
  
  void    operator()(MeasuredStar& m)
  {
    double val;
    bool   status = model_.Compute(m, pars_, val, deriv_);
    if(!status) return;
    
    double flux  = m.flux;
    double eflux = m.FluxSig();
    assert(eflux>0);
    
    double resid = sq( (flux-val)/eflux );
    chi2_ += resid;
    nmes_++;
    
    resids_.push_back(resid);
  }
  
  void    Finalize() 
  {
    if(!cmpsig_) return;
    unsigned int i, sz=resids_.size();
    double* residuals = new double[sz];
    for(i=0;i<sz;i++) residuals[i] = resids_[i];
    double mean=0;
    medres_ = sigres_ = 0.;
    Dmean_median_sigma(residuals, sz, mean, medres_, sigres_);
    delete[] residuals;
  }
  
  double  Chi2() const { return chi2_; }
  double  MedianRes() const { return medres_; }
  double  SigmaRes() const { return sigres_; }
  int     NMes() const { return nmes_; }
  
private:
  double chi2_;
  double medres_;
  double sigres_;
  
  Model  model_;
  Vect const&  pars_;
  Vect   deriv_;
  int nmes_;
  std::vector<double> resids_;
  bool cmpsig_;
};


template<class Model>
bool PhotFitter<Model>::FlipFlopStep()
{
  cout << "[PhotFitter<Model>::FlipFlopStep]" << endl;
  cleanChi2_ = false;
  
  unsigned int i;
  int igroup, ngroups = model_.NGroups();
  for(igroup=0;igroup<ngroups;igroup++)
    {
      unsigned int sz = model_.NPar(igroup);
      Mat  M(sz,sz);
      Vect B(sz);
      cout << "   *** igroup=" << igroup 
       << ", npar=" << sz << endl;
      
      MatrixFill<Model> fill(pars_, M, B, model_, igroup);
      LoopOnMeasurements(fill); 
      
      //      M.writeFits("toto.fits");
      //      exit(0);
      
      int status = cholesky_solve(M, B, "U");
      if(status!=0)
    {
      cout << "[PhotFitter<>::InitFit] cholesky_solve ERROR, status=" 
           << status << endl;
      M.writeFits("toto.fits");
      return false;
    }
      //      else 
      //    {
      //      cout << "[PhotFitter<>::InitFit] cholesky_solve OK, status="
      //           << status << endl;
      //    }
      
      // and now, update the --valid-- parameters
      //      unsigned int np = 0;
      unsigned int idx, ig = model_.GroupIndex(igroup);
      for(i=0;i<sz;i++)
    if( fill.IsValid(i) )
      {
        //      assert( !isinf(B(i)) ); // ceinture et bretelles
        idx = ig + i;
        pars_(idx) += B(i);
        //      cout << " par(" << idx << ") --> "
        //       << pars_(idx) << "(dpar=" << B(i) << ")" << endl;
        //      np++;
      }
        else 
      {
        cout << "Cannot update par(" << i << ")=" 
         << pars_(i) 
         << " B(" << i << ")=" << B(i)
         << endl;
      }
    }
  
  return true;
}


template<class Model>
bool PhotFitter<Model>::NewtonRaphsonStep()
{
  cleanChi2_ = false;
  
  unsigned int sz = model_.NPar();
  Mat  M(sz,sz);
  Vect B(sz);
  
  clock_t start;
  clock_t stop;
  
  start = clock();
  MatrixFill<Model> fill(pars_, M, B, model_, -1);
  LoopOnMeasurements(fill);
  stop = clock();
  cout << "[PhotFitter<>::NewtonRaphsonStep] " << (double)(stop-start) / (double)CLOCKS_PER_SEC
       << " to fill the fit matrix" << endl;
  
  start = clock();
  int status = cholesky_solve(M, B, "U");
  stop = clock();
  cout << "[PhotFitter<>::NewtonRaphsonStep] " << (double)(stop-start) / (double)CLOCKS_PER_SEC
       << " to invert it (using cholesky_solve)" << endl;
  
  if(status!=0)
    {
      cout << "[PhotFitter<>::NewtonRaphsonStep] cholesky_solve ERROR, status="
       << status
       << endl;
      M.writeFits("NewtonRaphson_DEBUG.fits");
      return false;
    }
  
  //  M.writeFits("NewtonRaphson_DEBUG.fits");
  
  // update the valid parameters
  unsigned int i;
  for(i=0;i<sz;i++)
    if( fill.IsValid(i) )
      {
    pars_(i) += B(i);
      }
    else
      {
    // should return false here ?????
    cout << "Cannot update pars(" << i << ")=" << pars_(i)
         << " <-> " << B(i)
         << endl;
      }

  return true;
}


template<class Model>
bool PhotFitter<Model>::Minimize(unsigned int maxiter)
{
  bool   status;
  double oldchi2;
  double chi2 = 99.E99;
  
  unsigned int i=0;
  while( 1 )
    {
      NewtonRaphsonStep();
      oldchi2 = chi2;
      chi2 = GetChi2();
      
      if(verbose_)
    {
      cout << "[PhotFitter<Model>::Minimize()] STEP #" << i
           << " CHI2="  << chi2;
    }
      
      if(oldchi2 < chi2)
    {
      // we should do something clever here. 
      // backtrack ?
      cout << endl
           << "[PhotFitter<>::Mininize()] ERROR increasing chi2 "
           << oldchi2 << " --> " << chi2 << endl;
      status = false;
      break;
    }
      
      if( (oldchi2-chi2) < minstep_ )
    {
      if(verbose_)
        {
          cout << " ... CONVERGED." << endl;
        }
      status = true;
      break;
    }
      
      if(i>=maxiter)
    {
      cout << endl
           << "[PhotFitter<>::Minimize()] ERROR, maxiter=" 
           << maxiter << " reached wo convergence." << endl;
      status = false;
      break;
    }
      
      cout << endl;
      i++;
    }
  
  status = ComputeCovMat();
  status = ComputeCorrMat();
  
  if(status)  updateFittedStars_();
  
  return status;
}


template<class Model>
bool PhotFitter<Model>::ComputeCovMat()
{
  epars_.allocate(model_.NPar());
  covmat_.allocate(model_.NPar(), model_.NPar());
  
  // we use epars temporarily
  MatrixFill<Model> fill(pars_, covmat_, epars_, model_, -1);
  LoopOnMeasurements(fill);
  
  int status = cholesky_solve(covmat_, epars_, "U");
  if(status!=0)
    {
      cout << "[PhotFitter<Model>::ComputeCovMat] cholesky_solve ERROR, status="
       << status << endl;
      covmat_.Zero();
      epars_.Zero();
      return false;
    }
  status = cholesky_invert(covmat_, "U");
  if(status!=0)
    {
      cout << "PB !";
      covmat_.Zero();
      epars_.Zero();
      return false;
    
    }
  
  unsigned int i;
  for(i=0;i<npar_;i++)
    epars_(i) = sqrt(covmat_(i,i));
  
  return true;
}


template<class Model>
bool PhotFitter<Model>::ComputeCorrMat()
{
  bool status = true;
  
  // this is a pretty lame test.
  // add an internal state record
  if( covmat_.SizeX() != model_.NPar() || 
      covmat_.SizeY() != model_.NPar() )
    status = ComputeCovMat();
  
  if(!status) return status;
  
  unsigned int i;
  unsigned int j;
  corrmat_.allocate(model_.NPar(), model_.NPar());
  
  for(i=0;i<model_.NPar();i++)
    {
      double cii = sqrt(covmat_(i,i));
      for(j=i;j<model_.NPar();j++)
    {
      double cjj = sqrt(covmat_(j,j));
      corrmat_(i,j) = covmat_(i,j) / (cii*cjj);
    }
    }
  
  return status;
}



template<class Model>
class OutlierRemover 
{
public:
  OutlierRemover(Model const& model, Vect const& pars, 
         double cut)
    : cut_(cut), nkilled_(0), model_(model), pars_(pars)
  {
    deriv_.allocate(model_.NPar());
  }
  
  ~OutlierRemover() { }
  
  void operator()(MeasuredStar& m)
  {
    double val;
    bool status = model_.Compute(m, pars_, val, deriv_);
    if(!status)
      {
    m.SetValid(false);
    nkilled_++;
      }
    
    double flux = m.flux;
    double eflux = m.FluxSig();
    assert(eflux>0);
    
    double resid = sq( (flux-val) / eflux );
    if(resid > cut_) 
      {
    m.SetValid(false);
    nkilled_++;
      }
  }
  
  void Finalize() { }
  
  unsigned int NKilled() const { return nkilled_; }
  
private:
  double cut_;
  unsigned int nkilled_;
  Model const& model_;
  Vect const& pars_;
  Vect deriv_;
};



template<class Model>
int PhotFitter<Model>::RemoveOutliers(double nsig)
{
  Chi2Acc<Model> acc(model_, pars_, true);
  LoopOnMeasurements(acc);
  double med = acc.MedianRes();
  double sig = acc.SigmaRes();
  double cut = med + nsig * sig;
  assert(cut>0.);
  
  cout << " med=" << med
       << " sig=" << sig
       << " cut=" << cut << endl;
  
  cout << "[PhotFitter<>::RemoveOutliers()] BEFORE OUTLIER REMOVAL: " << GetChi2() << endl;
  OutlierRemover<Model> rem(model_, pars_, cut);
  LoopOnMeasurements(rem);
  
  cleanChi2_ = false;
  cout << "[PhotFitter<>::RemoveOutliers()] " << rem.NKilled() << " measurements removed" << endl;
  cout << "[PhotFitter<>::RemoveOutliers()] AFTER OUTLIER REMOVAL: " << GetChi2() << endl;
  
  return rem.NKilled();
}



template<class Model>
unsigned int PhotFitter<Model>::GetNFittedPars() const
{
  unsigned int ret = 0;
  unsigned int i;
  for(i=0;i<model_.NPar();i++)
    //    if( !isinf(pars_(i)) )
    if( pars_(i)>=0 )
      ret++;
  
  return ret;
}

template<class Model>
double PhotFitter<Model>::GetChi2() const
{
  if(cleanChi2_) return chi2_;
  Chi2Acc<Model> acc(model_, pars_);
  LoopOnMeasurements(acc);
  chi2_ = acc.Chi2();
  nmes_ = acc.NMes();
  ndof_ = acc.NMes() - model_.NPar();
  cleanChi2_ = true;
  return chi2_;
}


template<class Model>
int PhotFitter<Model>::GetNMes() const
{
  if(cleanChi2_) return nmes_;
  GetChi2();
  return nmes_;
}

template<class Model>
int PhotFitter<Model>::GetNDoF() const
{
  if(cleanChi2_) return ndof_;
  GetChi2();
  return ndof_;
}


template<class Model>
Vect const& PhotFitter<Model>::GetPar() const
{
  return pars_;
}


template<class Model>
Vect const& PhotFitter<Model>::GetParError() const
{
  return epars_;
}


template<class Model>
Mat const& PhotFitter<Model>::GetCovMatrix() const
{
  return covmat_;
}

template<class Model>
Mat const& PhotFitter<Model>::GetCorrMatrix() const
{
  return corrmat_;
}


template<class Model>
void PhotFitter<Model>::updateFittedStars_()
{
  FittedStarIterator I;
  for(I=fittedstarstd::list_.begin();I!=fittedstarstd::list_.end();I++)
    {
      FittedStar* fs = (*I);
      if(!fs) continue;
      model_.UpdateFittedStar(*fs, pars_, epars_);
    }
}


template<class Model>
void PhotFitter<Model>::dumpParNTuple_(std::string const& ntname) const
{
  unsigned int i;
  ofstream ofs(ntname.c_str());
  ofs << "# num : " << endl
      << "# val : " << endl
      << "# err : " << endl
      << "# end"    << endl;
  
  for(i=0;i<model_.NPar();i++)
    {
      ofs << i << " "
      << pars_(i) << " "
      << epars_(i) << endl;
    }
  ofs.close();
}



// CHECK: WHAT ABOUT THE BAND ????
template<class Model>
void PhotFitter<Model>::dumpFittedStarNTuple_(const std::string& ntname) const
{
  ofstream ofs(ntname.c_str());
  ofs << "# index : " << endl
      << "# ra : " << endl
      << "# dec : " << endl
      << "# flux : " << endl
      << "# eflux : " << endl
      << "# nmeas : " << endl
      << "# ref : " << endl
      << "#end" << endl;
  
  FittedStarCIterator I;
  for(I=fittedstarstd::list_.begin();I!=fittedstarstd::list_.end();I++)
    {
      const FittedStar* fs = (*I);
      if(!fs) continue;
      const RefStar* rs = fs->GetRefStar();
      ofs << fs->Index() << " "
      << fs->X() << " "
      << fs->Y() << " "
      << fs->Flux() << " "
      << fs->FluxErr() << " "
      << fs->MeasurementCount() << " "
      << (rs!=0) << " "
      << endl;
    }
  
  ofs.close();
}




class MeasurementDumper
{
public:
  MeasurementDumper(ofstream& ofs, FittedStarMap const& fsm,
            unsigned int nusrvals=0)
    : ofs_(ofs), fsm_(fsm), nusrvals_(nusrvals)
  {
    writeHeader_(ofs_);
  }
  
  ~MeasurementDumper() { }
  
  void operator()(MeasuredStar const& ms)
  {
    CcdImage const& ccdim = ms.GetCcdImage();
    FittedStar const* fs = ms.GetFittedStar();
    double fluxfit = -1;
    double efluxfit = -1;
    double xdeg = -1;
    double ydeg = -1;
    int nmeas = -1;
    int nccds = -1;
    if(fs)
      {
    fluxfit = fs->Flux();
    efluxfit = fs->FluxErr();
    xdeg = fs->x;
    ydeg = fs->y; // CHECK THIS -- no ra,dec, but CTP coords.
    nmeas = fs->MeasurementCount();
    nccds = fsm_.NCcds(fs);
      }
    
    ofs_ << ms.x << " "
     << ms.y << " "
     << ms.flux << " "
     << ms.eflux << " "
     << fluxfit << " "
     << efluxfit << " "
     << xdeg << " "
     << ydeg << " "
     << ccdim.Shoot() - 700000 << " "
     << ccdim.Chip() << " "
     << ccdim.BandIndex() << " "
     << ccdim.AirMass() << " "
     << nmeas << " "
     << nccds << " ";
    
    // now, possibly the Monte-Carlo Fluxes ... 
    BaseStarWithError const* mcstar = ms.GetMCMeas();
    if(mcstar)
      ofs_ << mcstar->Flux() << " "
       << mcstar->EFlux() << " ";
    else 
      ofs_ << 0.0 << " " 
       << 0.0 << " ";
    
    BaseStarWithError const* mcstar_truth = ms.GetMCTruth();
    if(mcstar_truth)
      ofs_ << mcstar_truth->Flux();
    else
      ofs_ << 0.0 << " ";
    
    // and now, the usrVals
    assert(nusrvals_ <= ms.usrVals.size());
    unsigned int i;
    for(i=0;i<nusrvals_;i++)
      ofs_ << ms.usrVals[i] << " ";
    
    ofs_ << endl;
    
    //    cout << " ++++   >>>> " 
    //   << ms.GetMCTruth() << " "
    //   << ms.GetMCMeas() 
    //   << endl;
  }
  
  void Finalize() {}
  
private:
  ofstream& ofs_;
  FittedStarMap const& fsm_;
  unsigned int nusrvals_;

  void  writeHeader_(ofstream& ofs)
  {
    ofs << "# x : x in ccd" << endl
    << "# y : y in ccd" << endl
    << "# flux : flux measured on ccd" << endl
    << "# eflux : flux error" << endl
    << "# fluxfit : fitted flux" << endl
    << "# efluxfit : error on fitted flux" << endl
    << "# xdeg : x in degrees in tangent plane" << endl
    << "# ydeg : y in degrees in tangent plane" << endl
    << "# image : image number - 700000" << endl
    << "# ccd : ccd" << endl
    << "# band : band used for this measurement" << endl
    << "# airmass : airmass of this measurement" << endl
    << "# nm : number of measurements" << endl
    << "# nccd : number of diff ccds" << endl
    << "# mflux : true (distorted) flux on CCD" << endl
    << "# meflux : true flux error on CCD" << endl
    << "# iflux : initial flux" << endl;
    unsigned int i;
    for(i=0;i<nusrvals_;i++)
      ofs << "# usr" << i << " : " << endl;
    ofs << "# end" << endl;
  }
};



template<class Model>
void PhotFitter<Model>::dumpMeasuredStarNTuple_(const std::string& ntname,
                        unsigned int nusrvals) const
{
  ofstream ofs(ntname.c_str());
  FittedStarMap fsm(imstd::list_);
  MeasurementDumper dump(ofs, fsm, nusrvals);
  LoopOnMeasurements(dump);
  ofs.close();
}



template<class Model>
void PhotFitter<Model>::DumpNTuples(const std::string& dirname) const
{
  if(verbose_)
    {
      cout << "[PhotFitter<>::DumpNTuple(" << dirname << ")]" << endl;
    }
  
  if( !FileExists(dirname) )
    if( !MKDir(dirname.c_str()) )
      {
    cout << "[PhotFitter<>::DumpNTuple()] ERROR "
         << " unable to create dir=" << dirname
         << endl;
    return;
    }
  
  // the fit parameters
  std::string parntname = dirname + "/fitpars.ntuple";
  dumpParNTuple_(parntname);
  
  std::string fstarntname = dirname + "/fstars.std::list";
  dumpFittedStarNTuple_(fstarntname);
  
  std::string mstarntname = dirname + "/mstars.std::list";
  dumpMeasuredStarNTuple_(mstarntname);
  
  return;
}

#endif