photomfit2.h

// -*- C++ -*-
#ifndef PHOTOMFIT2_H
#define PHOTOMFIT2_H

#include <map>
#include <assert.h>
#include <iostream>
#include <countedref.h>

#include "refzp.h"
#include "ccdimage.h"
#include "measuredstar.h"
#include "fittedstar.h"



class Associations;


class PhotomModel : public RefCount 
{
public:
  virtual ~PhotomModel() {}
  
  virtual bool             isRef(CcdImage const&) const=0;
  virtual unsigned int     NScaleFactors() const=0;
  virtual unsigned int     PhotFactIndex(CcdImage const&) const=0;
  virtual double           GetPhotFact(CcdImage const&) const=0;
  virtual std::vector<double>   GetPhotFactors() const = 0;
  virtual std::vector<double>   GetPhotFactErrors() const = 0;
  
  virtual double           GetPhotFact(unsigned int index) const=0;
  virtual void             SetPhotFact(unsigned int index, double fact, double efact)=0;
  virtual unsigned int     CcdNum(unsigned int index) const=0;
};


// class OneGlocalScalePerShoot
// class OneGlocalScalePerCcdAndShoot

class OneGlobalScalePerCCD : public PhotomModel
{
public:
  OneGlobalScalePerCCD(std::map<unsigned int,double>& refscales);
  virtual ~OneGlobalScalePerCCD() {}
  
  virtual bool             isRef(CcdImage const&) const;
  virtual unsigned int     NScaleFactors() const;
  virtual unsigned int     PhotFactIndex(CcdImage const&) const;
  virtual double           GetPhotFact(CcdImage const&) const;
  virtual std::vector<double>   GetPhotFactors() const;
  virtual std::vector<double>   GetPhotFactErrors() const;
  
  virtual double           GetPhotFact(unsigned int index) const;
  virtual void             SetPhotFact(unsigned int index, double fact, double efact);
  virtual unsigned int     CcdNum(unsigned int index) const;
  
private:
  void            init();
  unsigned int    nscfact_;
  std::vector<unsigned int> refccds_;
  std::map<unsigned int,double> refscales_;
  std::vector<double>  scalefactors_; // size = 36
  std::vector<double>  scalefacterrors_; // size = 36
  std::map<unsigned int, unsigned int> ccd_; // ccd[index]
  std::map<unsigned int, unsigned int> indexes_; // index[ccd]
};



class PhotomFit2 
{
public:
  
  PhotomFit2(Associations& A, PhotomModel& model);
  
  bool                  Minimize(const int MaxIter);
  
  void                  RemoveOutliers(const double NSigCut);
  
  CcdImageList const&   GetCcdImageList() const { return ccdImageList_; }
  
  FittedStarList const& GetFittedStarList() const;
  
  double                ComputeChi2();
  
  PhotomModel&          Model() { return model_; }
  PhotomModel const&          Model() const { return model_; }

  // Get RID OF THIS. USELESS. -- or debug it
  int                   NValidFittedStars() const { return nfstars_; }
  
  void                  DumpFluxScales(const std::string& pattern) const;
  
  void                  SaveResults(const std::string& pattern) const;
  
private:
  void                  init();
  
  void                  update_star_fluxes();
  
  void                  update_scale_factors();
  
  int                   used_fittedstars();
  
  template<class Op>
  void   loop_on_measurements(CcdImageList&, Op& /* , bool skip_uninitialized_stars=false */);
  
  
  Associations&    assoc_;
  PhotomModel&     model_;
  FittedStarList&  fittedStarList_;
  
  CcdImageList     refCcds_;
  CcdImageList     otherCcds_;
  CcdImageList     ccdImageList_;
  
  unsigned int     nfstars_;
  unsigned int     nfmeasurements_;
  unsigned int     nrefmeasurements_;
  unsigned int     nfluxscales_;
  
  double   refChi2_;
  double   measuredChi2_;
};


#endif