Commit 8f06d719 authored by Reza  ANSARI's avatar Reza ANSARI
Browse files

Suite codage/restructuration programmes d'ajustement des auto-correlations et...

Suite codage/restructuration programmes d'ajustement des auto-correlations et cross-correlations, Reza 17/02/2019
parent 63b47b0a
......@@ -10,70 +10,69 @@
#include "sunitpcst.h"
#include "cxbeam.h"
#include "trkfit.h"
class My6CxSignalsB {
public:
My6CxSignalsB(vector< vector<double> > & v_time_data_, vector< vector< vector< complex<double> > > > & vv_cxdata_,
vector< vector< vector<double> > > & vv_cxerr_, // vector<double> & v_freqs_,
vector< SLinInterp1D > & v_interp_elev_, vector< SLinInterp1D > & v_interp_sazim_,
vector< vector<CxBeam> > & v_beams_, vector< vector<double> > & v_A_) // vector< vector< complex<double> > > & v_B_
: n_cxcor(6), v_time_data(v_time_data_), vv_cxdata(vv_cxdata_), vv_cxerr(vv_cxerr_), // v_freqs(v_freqs_),
v_interp_elev(v_interp_elev_), v_interp_sazim(v_interp_sazim_),
v_beams(v_beams_), v_A(v_A_)
My6CxSignalsB(vector< AcxDataSet> & acxd, vector< TrackSet > & trks)
: n_cxcor(6), acxd_(acxd), trk_(trks)
{
if (acxd_.size() != trk_.size())
throw ParmError("My6CxSignalsB::My6CxSignalsB(acxd,trks)/ERROR acxd_.size() != trk_.size()");
//DBG cout << " My6CxSignalsB() "<< endl;
}
// get the j-th signal for cross-correlation kcx (0<=kcx<=5)
virtual int getDataSignal(size_t j, size_t kcx, vector< complex<double> > & sig)
// get the the signal in dataset i , corresponding to the j-th track for cross-correlation kcx (0<=kcx<=5)
virtual int getDataSignal(size_t i, size_t j, size_t kcx, vector< complex<double> > & sig)
{
//DEL sig.resize(v_time_data[j].size());
vector< complex<double> > & vcxdata = vv_cxdata[j][kcx];
vector< complex<double> > & vcxdata = acxd_[i].vv_cxdata[j][kcx];
sig = vcxdata;
//DEL for(size_t k=0; k<sig.size(); k++) sig[k]=vcxdata[k];
return (int)sig.size();
}
// get the j-th signal for cross-correlation kcx (0<=kcx<=5)
virtual TVector< complex<double> > getDataSignal(size_t j, size_t kcx)
virtual TVector< complex<double> > getDataSignal(size_t i, size_t j, size_t kcx)
{
vector< complex<double> >sig;
getDataSignal(j, kcx, sig);
getDataSignal(i, j, kcx, sig);
TVector< complex<double> > rvs(sig);
//DEL for(size_t k=0; k<sig.size(); k++) rvs(k)=sig[k];
return rvs;
}
// get the beam for the j-th signal and the cross-correlation kcx (0<=kcx<=5)
inline CxBeam & getBeam(size_t j, size_t kcx)
// get the beam for the i-th data set and the cross-correlation kcx (0<=kcx<=5)
inline CxBeam & getBeam(size_t i, size_t kcx)
{
vector<CxBeam> & beams = v_beams[j];
return beams[kcx];
return acxd_[i].v_cxbeams[kcx];
}
// get the beam for the j-th signal and the cross-correlation kcx (0<=kcx<=5)
virtual int getExpectedSignal(size_t j, size_t kcx, vector< complex<double> > & sig, double phase, Vector3d & baseline_shift)
// get the expected signal for trackset i , corresponding to the j-th track for cross-correlation kcx (0<=kcx<=5)
virtual int getExpectedSignal(size_t i, size_t j, size_t kcx, vector< complex<double> > & sig, double phase, Vector3d & baseline_shift)
{
//DBG cout<<"*DBG*A*My6CxSignalsB::getExpectedSignal() j="<<j<<" size="<<v_time_data[j].size()<<" kcx="<<kcx<<endl;
sig.resize(v_time_data[j].size());
vector<double> & vtime = v_time_data[j];
vector<double> & vtime = acxd_[i].v_time_data[j];
sig.resize(vtime.size());
// vector< complex<double> > & vcxdata = vv_cxdata[j][nac_];
// vector<double> & vcxerr = vv_cxerr[j][nac_];
vector<double> & lesA = v_A[j];
double A = lesA[kcx];
CxBeam beam = getBeam(j, kcx);
double A = acxd_[i].v_Acx[j][kcx];
complex<double> B = acxd_[i].v_Bcx[j][kcx];
CxBeam beam = getBeam(i, kcx);
beam.ShiftBaseline(baseline_shift);
for(size_t k=0; k<vtime.size(); k++) {
double elev=v_interp_elev[j].YInterp(vtime[k]);
for(size_t l=0; l<vtime.size(); l++) {
double elev=trk_[i].v_interp_elev[j].YInterp(vtime[l]);
double thetasrc=Angle(90.-elev,Angle::Degree).ToRadian();
double sazim=v_interp_sazim[j].YInterp(vtime[k]);
double sazim=trk_[i].v_interp_sazim[j].YInterp(vtime[l]);
while (sazim > 360.) sazim -= 360.;
double phisrcdeg = 90.-sazim;
if (phisrcdeg < 0.) phisrcdeg+=360.;
complex<double> cxb=beam.Value(thetasrc, Angle(phisrcdeg,Angle::Degree).ToRadian());
// on fait conjugue complexe - car les donnees PAON4 sont Vi conj(Vj) , ici, on calcule conj(Vi) * Vj
sig[k] = complex<double>(A*cos(phase),A*sin(phase))*conj(cxb);
sig[l] = complex<double>(A*cos(phase),A*sin(phase))*conj(cxb);
sig[l] += B;
//DBG cout << " *DBG* "<<k<<" Theta,Phi-src = " << Angle(thetasrc).ToDegree() << " , " << Angle(phisrc).ToDegree()
//DBG << " elev="<<elev<<" phisrc="<<Angle(phisrc).ToDegree()<<" --> acb= "<<acb<< " sig="<<sig[k]<<endl;
......@@ -114,23 +113,22 @@ public:
}
return;
}
// get the beam for the j-th signal and the cross-correlation kcx (0<=kcx<=5)
virtual int getExpectedSignal(size_t j, size_t kcx, vector< complex<double> > & sig, const double* parm)
// get the expected signal for trackset i , corresponding to the j-th track for cross-correlation kcx (0<=kcx<=5)
virtual int getExpectedSignal(size_t i, size_t j, size_t kcx, vector< complex<double> > & sig, const double* parm)
{
double phase;
Vector3d baseline_shift;
getFromFitParam(kcx, parm, phase, baseline_shift);
double A = parm[1+3*j];
complex<double> B = complex<double>(parm[2+3*j], parm[3+3*j]);
return getExpectedSignal(j, kcx, sig, phase, baseline_shift);
return getExpectedSignal(i, j, kcx, sig, phase, baseline_shift);
}
virtual TVector< complex<double> > getExpectedSignal(size_t j, size_t kcx, double phase, Vector3d & baseline_shift)
// get the expected signal for trackset i , corresponding to the j-th track for cross-correlation kcx (0<=kcx<=5)
virtual TVector< complex<double> > getExpectedSignal(size_t i, size_t j, size_t kcx, double phase, Vector3d & baseline_shift)
{
vector< complex<double> > sig;
//DBG cout << "*DBG*getExpectedSignal() : calling getExpectedSignal(j="<<j<<" ,sig,...)"<<endl;
//DBG debug_gacfit=1;
getExpectedSignal(j, kcx, sig, phase, baseline_shift);
getExpectedSignal(i, j, kcx, sig, phase, baseline_shift);
TVector< complex<double> > rvs(sig);
return rvs;
}
......@@ -139,17 +137,19 @@ public:
{
double xi2=0.;
npts=0;
for(size_t j=0; j<v_time_data.size(); j++) { // loop over trk data sets
for(size_t kcx=0; kcx<n_cxcor; kcx++) { // loop over the six cross-correlations
vector< complex<double> > sig;
npts += getExpectedSignal(j, kcx, sig, parm);
vector< complex<double> > & vcxdata = vv_cxdata[j][kcx];
vector<double> & verr = vv_cxerr[j][kcx];
for(size_t i=0; i<vcxdata.size(); i++) {
double xx=(sig[i].real()-vcxdata[i].real())/verr[i];
xi2 += (xx*xx);
xx=(sig[i].imag()-vcxdata[i].imag())/verr[i];
xi2 += (xx*xx);
for(size_t i=0; i<acxd_.size(); i++) {
for(size_t j=0; j<acxd_[i].v_time_data.size(); j++) { // loop over tracks in the data set
for(size_t kcx=0; kcx<n_cxcor; kcx++) { // loop over the six cross-correlations
vector< complex<double> > sig;
npts += getExpectedSignal(i, j, kcx, sig, parm);
vector< complex<double> > & vcxdata = acxd_[i].vv_cxdata[j][kcx];
vector<double> & verr = acxd_[i].vv_cxerr[j][kcx];
for(size_t l=0; l<vcxdata.size(); l++) {
double xx=(sig[l].real()-vcxdata[l].real())/verr[l];
xi2 += (xx*xx);
xx=(sig[l].imag()-vcxdata[l].imag())/verr[l];
xi2 += (xx*xx);
}
}
}
}
......@@ -158,26 +158,16 @@ public:
}
size_t n_cxcor; // Number of Cx correlations ( 6 for PAON4 )
vector< vector<double> > & v_time_data;
vector< vector< vector< complex<double> > > > & vv_cxdata;
vector< vector< vector<double> > > & vv_cxerr;
vector< SLinInterp1D > & v_interp_elev;
vector< SLinInterp1D > & v_interp_sazim;
vector< vector<CxBeam> > & v_beams;
vector< vector<double> > & v_A; // vecteur des ampltudes pour chaque cross-correlation - determine par le fit d'avant
// vector< vector< complex<double> > > & v_B; // vecteur des offsets pour chaque cross-correlation - determine par le fit d'avant
vector< AcxDataSet> & acxd_;
vector< TrackSet > & trk_;
};
//-----------------------------------------------------------------
class My6CxGenXi2B : public GeneralXi2, public My6CxSignalsB {
public:
My6CxGenXi2B(vector< vector<double> > & v_time_data_, vector< vector< vector< complex<double> > > > & vv_cxdata_,
vector< vector< vector<double> > > & vv_cxerr_, // vector<double> & v_freqs_,
vector< SLinInterp1D > & v_interp_elev_, vector< SLinInterp1D > & v_interp_sazim_,
vector< vector<CxBeam> > & v_beams_, vector< vector<double> > & v_A_)
My6CxGenXi2B(vector< AcxDataSet> & acxd, vector< TrackSet > & trks)
: GeneralXi2(12) , // 12 = 3 phases + 3*3 position shifts
My6CxSignalsB(v_time_data_, vv_cxdata_, vv_cxerr_, v_interp_elev_, v_interp_sazim_, v_beams_, v_A_)
My6CxSignalsB(acxd, trks)
{
}
......
......@@ -12,7 +12,7 @@ MYOLISTHERE = $(OBJ)/p4autils.o $(OBJ)/visip4reader.o $(OBJ)/visp4winreader.o
MYOLISTHEREFIT = $(OBJ)/trkfit.o
# Define our target list
all : rdvisip4 visi2ntac visi2dtacx visi2tmfreq tstutls p4conv2fits msvis2dt visiavg filt_blind ckconvphases vis2ra p4vdblist rdthermrfi visamm trkacfit tfm2dt
all : rdvisip4 visi2ntac visi2dtacx visi2tmfreq tstutls p4conv2fits msvis2dt visiavg filt_blind ckconvphases vis2ra p4vdblist rdthermrfi visamm trkacfit tfm2dt trkacxfit
clean :
rm -f $(EXE)/rdvisip4 $(EXE)/visi2ntac $(EXE)/visi2dtacx $(EXE)/visi2tmfreq $(EXE)/p4conv2fits $(EXE)/msvis2dt $(EXE)/visiavg \
......@@ -95,6 +95,15 @@ $(EXE)/trkacfit : $(OBJ)/trkacfit.o $(MYOLISTHERE) $(MYOLISTHEREFIT)
$(OBJ)/trkacfit.o : trkacfit.cc $(MYINCLISTHERE) $(MYINCLISTHEREFIT)
$(CXXCOMPILE) -o $(OBJ)/trkacfit.o trkacfit.cc
trkacxfit : $(EXE)/trkacxfit
echo '---trkacxfit made'
$(EXE)/trkacxfit : $(OBJ)/trkacxfit.o $(MYOLISTHERE) $(MYOLISTHEREFIT)
$(CXXLINK) -o $(EXE)/trkacxfit $(OBJ)/trkacxfit.o $(MYOLISTHERE) $(MYOLISTHEREFIT) $(SOPHYAEXTSLBLIST)
$(OBJ)/trkacxfit.o : trkacxfit.cc $(MYINCLISTHERE) $(MYINCLISTHEREFIT)
$(CXXCOMPILE) -o $(OBJ)/trkacxfit.o trkacxfit.cc
$(OBJ)/trkfit.o : trkfit.cc $(MYINCLISTHERE) $(MYINCLISTHEREFIT)
$(CXXCOMPILE) -o $(OBJ)/trkfit.o trkfit.cc
......
/*---
PAON4 analysis :
Determining antenna pointing and dish size fitting parameters ...
(source track corresponds to a satellite track (or planet or sun ...)
Does also fit of phases for the 6 cross-correlations
Uses gacfit.h and gcxfit.h acbeam.h cxbeam.h
R. Ansari , February 2019
Evolution of trkacfit.cc , intended to replace it
Example command to run :
csh> ../AnaPaon4/Objs/trkacxfit -D 4.5 -sdec 10 -out A21_fac.txt dt_A21,trk_A21_43057,1358,1410,1278.5 dt_A21,trk_A21_43055,1512,1555,1278.5
----------------------------------------------------------*/
#include <stdlib.h>
#include <string.h>
#include <iostream>
#include <fstream>
#include <cmath>
#include <complex>
#include "array.h"
#include "randr48.h"
#include "skymap.h"
#include "samba.h"
#include "fftwserver.h"
#include "strutilxx.h"
#include "ctimer.h"
#include "tarrinit.h"
#include "skymapinit.h"
#include "fitsioserver.h"
#include "slininterp.h"
#include "p4autils.h"
#include "acbeam.h"
#include "gacfit.h"
#include "gcxfit.h"
#include "gcxfitbaseline.h"
#include "trkfit.h"
using namespace std;
using namespace SOPHYA;
//---------------------------------------------------------------------------
//--- Grouping here beam & computation parameters -------------
//--- Global parameters
// static double latitude=45.; // latitude in degree
// static double longitude=0.; // longitude in degree
static double D_dish = 4.5; // effective dish size = effD*D
static bool fggaussbeam=true;
// --- if true, perform Beam normalisation
// static bool fgdobeamnormalise=false;
static string outfilename; // output filename for auto-correlation fit parameters
static string checkfilename; // output filename for PAON4 and computed auto-cor using fitted parameters
static bool fguseAac4Cx=true; //if true, use fitted Amplitudes on autocor for initial Cross-cor fit parameter value
// static bool fgoutfits = false; // true -> write output in fits format
static int prtlevel=0; // print level
static vector<string> trksetfiles; // datacard files defining the track/data sets
//--- End of list of global parameters
//-------- Declaration of functions in this file - code after the main, below
int decode_args(int narg, char* arg[]);
//--------------------------------------------------------------
//--------- Main program -------------
//--------------------------------------------------------------
int main (int narg, char* arg[])
{
cout << " ------ trkacxfit.cc : Fitting auto-correlations using track (satellites...) files ------- " << endl;
int rc = 0;
try {
SophyaInit();
int rcda=decode_args(narg, arg);
if (rcda) return rcda;
vector<AcxDataSet> v_acxd;
vector<TrackSet> v_trk;
for(size_t i=0; i<trksetfiles.size(); i++) {
TrkInputDataSet tids(trksetfiles[i]);
AcxDataSet acxd(tids);
TrackSet tks(tids);
ACxSetFitter acxfit(acxd, tks);
string acofile = "ac_"+outfilename;
acxfit.doACfit(acofile);
string acckfile = "ac_"+checkfilename;
acxfit.saveExpectedAC(acckfile);
string cxofile = "cx_"+outfilename;
string cxckfile = "ac_"+checkfilename;
acxfit.doCxfit(cxofile, cxckfile, fguseAac4Cx);
v_acxd.push_back(acxd);
v_trk.push_back(tks);
}
Timer tm("trkacxfit");
} // End of try bloc
catch (PThrowable & exc) {
cerr << " trkacxfit.cc: Catched Exception (PThrowable)" << (string)typeid(exc).name()
<< " - Msg= " << exc.Msg() << endl;
rc = 99;
}
catch (std::exception & e) {
cerr << " trkacxfit.cc: Catched std::exception "
<< " - what()= " << e.what() << endl;
rc = 98;
}
catch (...) {
; cerr << " trkacxfit.cc: some other exception (...) was caught ! " << endl;
rc = 97;
}
cout << " ------------ END OF trkacxfit.cc (Rc= " << rc << ") ------------ " << endl;
return rc;
}
//-----------------------------------------------------------------------------------------
//-----------------------------------------------------------------------------------------
/* --Fonction-- */
int decode_args(int narg, char* arg[])
//! Decode program arguments
{
bool fglonghelp=false;
if ((narg>1)&&(strcmp(arg[1],"-h")==0)) fglonghelp=true;
if ((narg<2)||fglonghelp) {
cout << " trkacxfit : fit array geometry and \n"
<< " Usage: trkacxfit [-options] track_Set_1 [track_Set_2] [track_Set_3 ...] \n"
<< " options: -D dish_diameter -ngb -out OutBaseFileName -ckf CheckBaseFileName \n"
<< " -prt PrintLevel \n"<<endl;
if (!fglonghelp) {
cout << " trkacxfit -h to get option description " << endl;
return 2;
}
cout << " -D dish_diameter : define effective dish diameter (in meter D*eff , def=4.5) \n"
<< " -ngb : Use non gaussian beam profile (Bessel j1(angle)) \n"
<< " -prt PrintLevel: specify print level \n"
<< " -out OutFileName : Output file (text) to save fitted parameters vectors (def trkfit.txt) \n"
<< " -ckf CheckFileName : Output file (PPF) to save expected signals (def chktrkfit.ppf) \n"<<endl;
return 1;
}
outfilename = "trkfit.txt";
checkfilename = "chktrkfit.ppf";
prtlevel=0;
D_dish = 4.5;
fggaussbeam=true;
vector<string> lastargs;
while (narg>1) {
string fbo = arg[1];
if (fbo=="-D") { // dish size (effective dish diameter)
if (narg<2) { cout << "trkacxfit/decode_args missing/bad argument, -h for help " << endl; return -1; }
D_dish=atof(arg[2]); arg+=2; narg-=2; lastargs.clear();
}
else if (fbo=="-ngb") { // Use Non gaussian beam
fggaussbeam=false; arg++; narg--; lastargs.clear();
}
else if (fbo=="-out") { // output file name
if (narg<2) { cout << "trkacxfit/decode_args missing/bad argument, -h for help " << endl; return -1; }
outfilename=arg[2]; arg+=2; narg-=2; lastargs.clear();
}
else if (fbo=="-ckf") { // output file name for checking (data & fitted
if (narg<2) { cout << "trkacxfit/decode_args missing/bad argument, -h for help " << endl; return -1; }
checkfilename=arg[2]; arg+=2; narg-=2; lastargs.clear();
}
else if (fbo=="-prt") { // print level
if (narg<2) { cout << "trkacxfit/decode_args missing/bad argument, -h for help " << endl; return -1; }
prtlevel=atoi(arg[2]); arg+=2; narg-=2; lastargs.clear();
}
else { arg++; narg--; lastargs.push_back(fbo); }
}
if (lastargs.size() < 1) {
cout << " trkacxfit/ERROR : No input track/data specified (trkacxfit -h for usage) \n";
return 9;
}
cout << " ------------------- trkacxfit/run parameters:"<<endl;
cout << " Beam: D_dish (initial guess)="<<D_dish<<" GaussianBeam ? "<<(fggaussbeam?"Yes":"No")<<endl;
cout << " OutFileName= "<<outfilename<<" CheckFileName= "<<checkfilename<<endl;
cout << " --- TrackSetFiles: (NbFiles="<<lastargs.size()<<endl;
trksetfiles = lastargs;
for (size_t i=0; i<trksetfiles.size(); i++) {
cout << "["<<i<<"] : "<< trksetfiles[i]<<endl;
}
cout<<"---------------------------------------------------------------------------------"<<endl;
return 0;
}
......@@ -331,15 +331,17 @@ size_t TrackSet::ReadData(TrkInputDataSet & tkds)
//------------------------ TrackSet -------------------------------------
ACxSetFitter::ACxSetFitter(AcxDataSet & data, TrackSet & tks)
: fggaussbeam_(true), D_dish(5.), acxd_(data), tks_(tks),
: fggaussbeam_(true), D_dish(5.), acxd_(data), tks_(tks), fit_ac_done(false), fit_cx_done(false),
v_RcFit_ac(tks.getNbAutoCor()), v_xi2red_ac(tks.getNbAutoCor()),
v_Ddish(tks.getNbAutoCor()), v_thetaant(tks.getNbAutoCor()),
v_phiant(tks.getNbAutoCor()), v_A(tks.getNbAutoCor()), v_B(tks.getNbAutoCor()),
v_err_Ddish(tks.getNbAutoCor()), v_err_thetaant(tks.getNbAutoCor()),
v_err_phiant(tks.getNbAutoCor()), v_err_A(tks.getNbAutoCor()), v_err_B(tks.getNbAutoCor()),
v_acbeams(tks.getNbAutoCor()),
v_RcFit_cx(tks.getNbCrossCor()), v_xi2red_cx(tks.getNbCrossCor()),
v_phase(tks.getNbCrossCor()), v_Acx(tks.getNbCrossCor()),
v_err_phase(tks.getNbCrossCor()), v_err_Acx(tks.getNbCrossCor())
v_err_phase(tks.getNbCrossCor()), v_err_Acx(tks.getNbCrossCor()),
v_cxbeams(tks.getNbCrossCor())
{
}
......@@ -446,6 +448,14 @@ int ACxSetFitter::doACfit(string outfilename)
ofr <<setw(8)<<A<<" "<<setw(8)<<err_A<<" "<<setw(8)<<B<<" "<<setw(8)<<err_B<<" ";
}
ofr << endl;
double clight = PhysQty::c().SIValue();
double lambda = clight/(acxd_.v_freqs[0]*1.e6);
ACBeam acb1(Dfit, thetaant, phiant, lambda);
acb1.setGaussianLobe(fggaussbeam_);
ACBeam acb2(Dfit, thetaant, phiant, lambda);
acb2.setGaussianLobe(fggaussbeam_);
Vector3d baseline0(0.,0.,0.);
v_acbeams[ii]=CxBeam(acb1, acb2, baseline0);
}
else {
cout << "---Fit failed for "<<ii+1<<"--- Fit_Error, rc = " << rcfit << " nstep="<<mFit.GetNStep()<<endl;
......@@ -454,7 +464,8 @@ int ACxSetFitter::doACfit(string outfilename)
}
}
fit_ac_done=true;
acxd_.v_acbeams=v_acbeams;
return 0;
}
......@@ -517,6 +528,10 @@ int ACxSetFitter::doCxfit(string outfilenamecx, string outcheckfilenamecx, bool
size_t Anum2[6]={1,2,3,2,3,3};
for(size_t ii=0; ii<NB_CXCORS; ii++) {
v_Acx[ii].resize(NTRK);
v_Bcx[ii].resize(NTRK);
for(size_t j=0; j<NTRK; j++) {
v_Acx[ii][j]=1.; v_Bcx[ii][j]=complex<double>(0.,0.);
}
Vector3d baseline=P4Coords::getBaseline(Anum1[ii]+1,Anum2[ii]+1);
cout << "--------- 1."<<ii+1<<" doCxfit() Doing fit for CrossCor= " << ii << " FxF= "
<< Anum1[ii]+1<<"x"<<Anum2[ii]+1<<" Baseline="<<baseline<<endl;
......@@ -536,6 +551,7 @@ int ACxSetFitter::doCxfit(string outfilenamecx, string outcheckfilenamecx, bool
ACBeam acb2(v_Ddish[Anum2[ii]], v_thetaant[Anum2[ii]], v_phiant[Anum2[ii]], lambda);
acb2.setGaussianLobe(fggaussbeam_);
CxBeam cxbeam(acb1, acb2, baseline);
v_cxbeams[ii]=cxbeam;
MyCxGenXi2 gxi2( acxd_.v_time_data, acxd_.vv_cxdata, acxd_.vv_cxerr,
tks_.v_interp_elev, tks_.v_interp_sazim, cxbeam, ii);
......@@ -639,7 +655,8 @@ int ACxSetFitter::doCxfit(string outfilenamecx, string outcheckfilenamecx, bool
double Ai=(useAac?Aac:v_amp[j]);
double A=mFit.GetParm(1+3*j); double err_A=mFit.GetParmErr(1+3*j);
cout << " Trk["<<j<<"] A= "<<A<<" +/- "<<err_A<<" (A/Ai="<<A/Ai<<")"<<endl;
v_Acx[ii][j]=A; // v_B[ii][j]=B;
v_Acx[ii][j]=A;
v_Bcx[ii][j]=complex<double>(0.,0.);
v_err_Acx[ii][j]=err_A;
ofr <<setw(8)<<A<<" "<<setw(8)<<err_A<<" ";
}
......@@ -659,5 +676,9 @@ int ACxSetFitter::doCxfit(string outfilenamecx, string outcheckfilenamecx, bool
}
}
if (pox) delete pox;
fit_cx_done=true;
acxd_.v_cxbeams=v_cxbeams;
acxd_.v_Acx=v_Acx;
acxd_.v_Bcx=v_Bcx;
return 0;
}
......@@ -15,6 +15,8 @@
#include "skymap.h"
#include "histats.h"
#include "cxbeam.h"
using namespace std;
using namespace SOPHYA;
......@@ -80,6 +82,10 @@ public:
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
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
};
......@@ -130,6 +136,8 @@ public:
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;
......@@ -138,12 +146,16 @@ public:
vector<double> v_err_Ddish;
vector<double> v_err_thetaant, v_err_phiant;
vector< vector<double> > v_err_A, v_err_B;
vector<double> v_phase;
vector< vector<double> > v_Acx;
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< vector<double> > v_Acx;
vector< vector< complex<double> > > v_Bcx; // Offset for the six cross-cors after fit
vector<double> v_err_phase;
vector< vector<double> > v_err_Acx;
vector< vector< complex<double> > > v_err_Bcx;
vector< CxBeam > v_cxbeams; // the six cross correlations after fit
};
//-----------------------------------------------------------------------
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment