trkfit.h 8.32 KB
Newer Older
1 2 3 4 5
/*  PAON4 analysis software 
    classes and functions to read in and perform array geometry determination 
    using satellites and celestial sources tracks  
    R. Ansari, Fevrier 2019                                             */

6 7 8
#ifndef TRKFIT_H_SEEN
#define TRKFIT_H_SEEN

9 10 11
#include <string>
#include <vector>
#include <iostream>
12 13 14 15 16
#include <complex>

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

18 19
#include "cxbeam.h"

20 21 22
using namespace std;
using namespace SOPHYA;

23 24
void TrkFit_SetPrintLevel(int lev=0);

25 26 27 28 29 30
//-------------------------------------------------------------------------
//------------------- TrkInputDataSet -------------------------------------
//-----   Class describing the input data set list  
class TrkInputDataSet {
public:
  TrkInputDataSet() : zenang(0.) , theta_0(0.) , phi_0(0.) { }
31 32 33
  TrkInputDataSet(string dcfilename, string inp_path="");    
  void setInputBasePath(string inp_path);
  // create the object reading filenames and other parameters from the datacard file
34 35 36 37 38 39 40 41 42 43 44
  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;
45 46
  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 
47 48 49 50 51 52
  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);  }
53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84

//-------------------------------------------------------------------------
//------------------------ 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;
85 86 87
  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    
88 89 90
  double zenang;     // zenith angle = shift with declination, used for initial value of fitted angles
  double theta_0, phi_0;             // Theta, phi angles corresponding

91
  vector< CxBeam > v_acbeams;     // the four aut-correlations beams after fit 
92
  vector< CxBeam > v_cxbeams;  // the six cross correlations after fit 
93 94 95
  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 
96 97
  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 
98 99 100 101 102 103 104 105 106
};


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

class TrackSet {
public:
107 108 109
  // 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);
  
110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139
  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  
140
  int doCxfit(string outfilenamecx, bool useAac=true, bool fgphi0only=true);  // outfilename : fit results   
141 142 143

  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  
144 145 146 147 148 149 150 151

  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_;
152 153
  bool fit_ac_done;  // true -> doACfit() called 
  bool fit_cx_done;  // true -> doACfit() called 
154 155 156 157 158 159 160 161
  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; 
162
  vector< CxBeam > v_acbeams;     // the four aut-correlations beams after fit 
163 164
  vector<int> v_RcFit_cx;
  vector<double> v_xi2red_cx;
165
  vector<double> v_phase;
166 167
  vector<double> v_phi_0;
  vector<double> v_a_phi;
168 169
  vector< vector<double> > v_Acx;
  vector< vector< complex<double> > > v_Bcx;   // Offset for the six cross-cors after fit 
170 171
  vector<double> v_err_phi_0;
  vector<double> v_err_a_phi;
172
  vector< vector<double> > v_err_Acx;
173 174
  vector< vector< complex<double> > > v_err_Bcx;   
  vector< CxBeam > v_cxbeams;  // the six cross correlations after fit 
175 176
};

177 178 179 180 181 182 183 184 185

//-----------------------------------------------------------------------
//-------------------------- 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);
186 187
  ~CxBaselineFitter();

188 189
  void initFitParams();

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

193
  int doCheck();
194 195
  int doSimplexMinimize();

196
  int saveExpectedCx(string outcheckfilename);  // outfilename PPF (or fits ?) file with the expected signals  
197 198 199 200 201 202

  vector<AcxDataSet> & v_acxd;
  vector<TrackSet> & v_trks;
  size_t tot_ntrks;
  
  bool fit_done;
203
  bool simplex_done;
204 205
  int rcfit;
  double xi2red;
206 207 208 209
  vector<double> v_phi_0;
  vector<double> v_err_phi_0;
  vector<double> v_a_phi;
  vector<double> v_err_a_phi;
210
  vector< Vector3d > v_baselineshits;
211 212 213
  vector< Vector3d > v_err_baselineshits;
  double * bestfitparam;
  double * err_bestfitparam;
214 215
};

216 217
//-----------------------------------------------------------------------
#endif