trkfit.h

/*  PAON4 analysis software 
    classes and functions to read in and perform array geometry determination 
    using satellites and celestial sources tracks  
    R. Ansari, Fevrier 2019                                             */

#ifndef TRKFIT_H_SEEN
#define TRKFIT_H_SEEN

#include <string>
#include <vector>
#include <iostream>
#include <complex>

#include "slininterp.h"
#include "skymap.h"
#include "histats.h"

#include "cxbeam.h"

using namespace std;
using namespace SOPHYA;

void TrkFit_SetPrintLevel(int lev=0);

//-------------------------------------------------------------------------
//------------------- TrkInputDataSet -------------------------------------
//-----   Class describing the input data set list  
class TrkInputDataSet {
public:
  TrkInputDataSet() : zenang(0.) , theta_0(0.) , phi_0(0.) { }
  TrkInputDataSet(string dcfilename, string inp_path="");    
  void setInputBasePath(string inp_path);
  // create the object reading filenames and other parameters from the datacard file
  size_t ReadDatacardFile(string dcfilename);
  ostream & Print(ostream & os) const;
  
  inline size_t  NbTrk() const { return dataflnm.size(); }
  
  string dcfilename_;
  string input_base_path;
  vector<string> dataflnm;    // PAON4 data Auto+Cross corr DataTable filenames 
  vector<double> tstart, tend;   // time interval in seconds (start,end time)  
  vector<double> v_freqs;     // Central frequency associated to each visibilty = f(time) dataset
  vector<string> trkflnm;
  vector<bool> v_noAC;     // if true -> dont use this track and associated visibilities for AutoCorrelation fit 
  vector<bool> v_noCx;     // if true -> dont use this track and associated visibilities for Cross-Correlation fit 
  double zenang;     // zenith angle = shift with declination, used for initial value of fitted angles
  double theta_0, phi_0;             // Theta, phi angles corresponding
};

inline ostream& operator << (ostream& s, TrkInputDataSet const & trk)
  {  trk.Print(s);  return(s);  }

//-------------------------------------------------------------------------
//------------------------ ACxDataSet -------------------------------------
//---  Class containing the Auto and Cross-correlation data for a track set 
 
class AcxDataSet {
public:
  AcxDataSet() : tot_npoints(0),zenang(0.),theta_0(0.),phi_0(0.) { }
  AcxDataSet(AcxDataSet const & a);
  AcxDataSet(TrkInputDataSet & tkds);
  
  size_t ReadData(TrkInputDataSet & tkds);

  AcxDataSet & operator = (AcxDataSet const & a);

  inline size_t getTotNPoints() const { return tot_npoints; }

  inline size_t  NbTrk() const { return v_time_data.size(); }

  static inline size_t getNbAutoCor() { return 4; }
  static inline size_t getNbCrossCor() { return 6; }

  vector< vector<double> > v_time_data;
  vector< vector< vector<double> > > vv_data;   
  vector< vector< vector<double> > > vv_err;   
  vector< vector<double> > v_min_data;
  vector< vector<double> > v_max_data;
  vector< vector< vector< complex<double> > > > vv_cxdata;   
  vector< vector< vector<double> > > vv_cxerr;   
  vector< vector<double> > v_min_cxdata;
  vector< vector<double> > v_max_cxdata;
  size_t tot_npoints;
  vector<double> v_freqs; 
  vector<bool> v_noAC;     // if true -> dont use this track and associated visibilities for AutoCorrelation fit 
  vector<bool> v_noCx;     // if true -> dont use this track and associated visibilities for Cross-Correlation fit    
  double zenang;     // zenith angle = shift with declination, used for initial value of fitted angles
  double theta_0, phi_0;             // Theta, phi angles corresponding

  vector< CxBeam > v_acbeams;     // the four aut-correlations beams after fit 
  vector< CxBeam > v_cxbeams;  // the six cross correlations after fit 
  vector<double> v_phase;    // fitted phases for individual cross-correlations at reference frequency
  vector<double> v_phi_0;    // fitted phases (Phi0) for individual cross-correlations 
  vector<double> v_a_phi;    // fitted phases (a_phi = slope with frequency) for individual cross-correlations 
  vector< vector<double> > v_Acx;   // Amplitudes for the six cross-cors after fit 
  vector< vector< complex<double> > > v_Bcx;   // Offset for the six cross-cors after fit 
};


//-----------------------------------------------------------------------
//-------------------------- TrackSet -----------------------------------
//-----      Class containing a set of tracks 

class TrackSet {
public:
  // Read a track file and return a vector of time (in second), vector of elevations, azimuth and the corresponding linear interpolation
  static size_t ReadTrackFile(string flnm, vector<double> & tims, vector<double> & elevs, vector<double> & azims, SLinInterp1D & li_elev, SLinInterp1D & li_sazim);
  
  TrackSet()  { }
  TrackSet(TrackSet const & a);
  TrackSet(TrkInputDataSet & tkds);
  
  size_t ReadData(TrkInputDataSet & tkds);

  TrackSet & operator = (TrackSet const & a);

  inline size_t  NbTrk() const { return v_time_sat.size(); }

  static inline size_t getNbAutoCor() { return 4; }
  static inline size_t getNbCrossCor() { return 6; }
  
  vector< vector<double> > v_time_sat;
  vector< vector<double>  > v_sat_elev;
  vector< vector<double>  > v_sat_azim;
  
  vector< SLinInterp1D > v_interp_elev;
  vector< SLinInterp1D > v_interp_sazim;
};  

//-----------------------------------------------------------------------
//-------------------------- ACxSetFitter -------------------------------
//--       Class for fitting a poitings (AutoCorrelation) and phases 
//--       on cross-correlations - 6 phases, each of the 6 cross-corr independtly 

class ACxSetFitter {
public:
  ACxSetFitter(AcxDataSet & data, TrackSet & tks);
  int doACfit(string outfilename);  // outfilename : fit results  
  int doCxfit(string outfilenamecx, bool useAac=true, bool fgphi0only=true);  // outfilename : fit results   

  int saveExpectedAC(string outcheckfilename);  // outfilename PPF (or fits ?) file with the expected AC signals  
  int saveExpectedCx(string outcheckfilename);  // outfilename PPF (or fits ?) file with the expected signals  

  inline void setD_dish(double D=5.) { D_dish=D; }
  inline void setGaussianBeam(bool fg=true) { fggaussbeam_=fg; }

  bool fggaussbeam_;
  double D_dish;
  AcxDataSet & acxd_;
  TrackSet & tks_;
  bool fit_ac_done;  // true -> doACfit() called 
  bool fit_cx_done;  // true -> doACfit() called 
  vector<int> v_RcFit_ac;
  vector<double> v_xi2red_ac;
  vector<double> v_Ddish;
  vector<double> v_thetaant, v_phiant;
  vector< vector<double> > v_A, v_B; 
  vector<double> v_err_Ddish;
  vector<double> v_err_thetaant, v_err_phiant;
  vector< vector<double> > v_err_A, v_err_B; 
  vector< CxBeam > v_acbeams;     // the four aut-correlations beams after fit 
  vector<int> v_RcFit_cx;
  vector<double> v_xi2red_cx;
  vector<double> v_phase;
  vector<double> v_phi_0;
  vector<double> v_a_phi;
  vector< vector<double> > v_Acx;
  vector< vector< complex<double> > > v_Bcx;   // Offset for the six cross-cors after fit 
  vector<double> v_err_phi_0;
  vector<double> v_err_a_phi;
  vector< vector<double> > v_err_Acx;
  vector< vector< complex<double> > > v_err_Bcx;   
  vector< CxBeam > v_cxbeams;  // the six cross correlations after fit 
};


//-----------------------------------------------------------------------
//-------------------------- CxBaselineFitter -------------------------------
//--       Class for fitting baselines and phases using the 6 cross correlations  
//--       on a set of tracks 

class CxBaselineFitter {
public:
  CxBaselineFitter(vector<AcxDataSet> & v_data, vector<TrackSet> & v_tks);
  ~CxBaselineFitter();

  void initFitParams();

// if  fgfixbaseline = true, fit phases only  , if fgphi0only : phases independent of frequency , aphi=0
  int dofit(string outfilename, bool fgfixbaseline=false, bool fgphi0only=true);  

  int doCheck();
  int doSimplexMinimize();

  int saveExpectedCx(string outcheckfilename);  // outfilename PPF (or fits ?) file with the expected signals  

  vector<AcxDataSet> & v_acxd;
  vector<TrackSet> & v_trks;
  size_t tot_ntrks;
  
  bool fit_done;
  bool simplex_done;
  int rcfit;
  double xi2red;
  vector<double> v_phi_0;
  vector<double> v_err_phi_0;
  vector<double> v_a_phi;
  vector<double> v_err_a_phi;
  vector< Vector3d > v_baselineshits;
  vector< Vector3d > v_err_baselineshits;
  double * bestfitparam;
  double * err_bestfitparam;
};

//-----------------------------------------------------------------------
#endif