Commit 080707b4 authored by Reza  ANSARI's avatar Reza ANSARI
Browse files

Suite implementation/debug fit baselines, Reza 19/02/2019

parent 3517af9b
......@@ -31,6 +31,13 @@ public:
v_lambdas[k] = clight/(v_freqs[k]*1.e6);
}
virtual TVector<double> getTimeVec(size_t j)
{
vector<double> & tvec = v_time_data[j];
TVector<double> tmv(tvec);
return tmv;
}
virtual int getDataSignal(size_t j, vector<double> & sig)
{
//DEL sig.resize(v_time_data[j].size());
......
......@@ -24,6 +24,13 @@ public:
//DBG cout << " MyCxSignal() CxCorr=ncx "<<ncx_<< endl;
}
virtual TVector<double> getTimeVec(size_t j)
{
vector<double> & tvec = v_time_data[j];
TVector<double> tmv(tvec);
return tmv;
}
virtual int getDataSignal(size_t j, vector< complex<double> > & sig)
{
//DEL sig.resize(v_time_data[j].size());
......
......@@ -22,6 +22,14 @@ public:
//DBG cout << " My6CxSignalsB() "<< endl;
}
// get the the TimeVec in dataset i , corresponding to the j-th track
virtual TVector<double> getTimeVec(size_t i, size_t j)
{
vector<double> & tvec = acxd_[i].v_time_data[j];
TVector<double> tmv(tvec);
return tmv;
}
// 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)
{
......@@ -133,6 +141,17 @@ public:
return rvs;
}
// 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, const double* parm)
{
vector< complex<double> > sig;
//DBG cout << "*DBG*getExpectedSignal() : calling getExpectedSignal(j="<<j<<" ,sig,...)"<<endl;
//DBG debug_gacfit=1;
getExpectedSignal(i, j, kcx, sig, parm);
TVector< complex<double> > rvs(sig);
return rvs;
}
virtual double getXi2(const double* parm, int & npts)
{
double xi2=0.;
......
......@@ -129,8 +129,9 @@ int main (int narg, char* arg[])
sprintf(fnbuff,"cx_%d_",(int)(i+1));
string cxofile = fnbuff+outfilename;
string cxckfile = fnbuff+checkfilename;
acxfit.doCxfit(cxofile, cxckfile, fguseAac4Cx);
acxfit.doCxfit(cxofile, fguseAac4Cx);
acxfit.saveExpectedCx(cxckfile);
v_acxd.push_back(acxd);
v_trk.push_back(tks);
}
......@@ -139,7 +140,8 @@ int main (int narg, char* arg[])
CxBaselineFitter cxbfit(v_acxd, v_trk);
string cxbofile = "cxb6_"+outfilename;
string cxbckfile = "cxb6_"+checkfilename;
cxbfit.dofit(cxbofile, cxbckfile);
cxbfit.dofit(cxbofile);
cxbfit.saveExpectedCx(cxbckfile);
}
Timer tm("trkacxfit");
} // End of try bloc
......
......@@ -520,13 +520,18 @@ int ACxSetFitter::saveExpectedAC(string outcheckfilename)
Vector expsignal = macs.getExpectedSignal(j, Ddishfit, thetafit, phifit, A, B);
sprintf(oname,"simac_%d_%d",(int)ii+1,(int)j+1);
pos << PPFNameTag(oname)<<expsignal;
if (ii==0) {
Vector tmvec = macs.getTimeVec(j);
sprintf(oname,"tim_%d",(int)j+1);
pos << PPFNameTag(oname)<<tmvec;
}
}
}
return 0;
}
int ACxSetFitter::doCxfit(string outfilenamecx, string outcheckfilenamecx, bool useAac)
int ACxSetFitter::doCxfit(string outfilenamecx, bool useAac)
{
size_t NB_ANTENNES=acxd_.getNbAutoCor(); // nombre d'antennes
size_t NB_CXCORS=acxd_.getNbCrossCor();
......@@ -535,11 +540,11 @@ int ACxSetFitter::doCxfit(string outfilenamecx, string outcheckfilenamecx, bool
cout << "======================================================================================"<<endl;
cout << "---------- ACxSetFitter::doCxfit() ; Performing cross-cor phase fit ..."<<endl;
if (useAac) cout << " ... Using Amplitude from auto-correlations fit for initial fit parameter value..."<<endl;
POutPersist * pox = NULL ;
/* DEL POutPersist * pox = NULL ;
if (outcheckfilenamecx.length()>0) {
cout << "... expected cross-cor for fitted params (and cadata) will be saved to file "<<outcheckfilenamecx<<endl;
pox = new POutPersist(outcheckfilenamecx);
}
} */
ofstream ofr(outfilenamecx.c_str());
ofr << "#### cross-cor phase fit (ACxSetFitter::doCxfit() ) "<<endl
<< "## NumCxCor Xi2red Phase err_Phase (deg) A0 err_A0 A1 err_A1 ..."<<endl;
......@@ -593,10 +598,10 @@ int ACxSetFitter::doCxfit(string outfilenamecx, string outcheckfilenamecx, bool
Vector asig = SOPHYA::abs(signal);
double mins, maxs;
asig.MinMax(mins, maxs);
if (pox) {
/* if (pox) {
sprintf(oname,"cx_%d_%d",(int)ii+1,(int)j+1);
(*pox) << PPFNameTag(oname)<<signal;
}
} */
TVector< complex<double> > expsignal = gxi2.getExpectedSignal(j, 0., A);
Vector aexpsig = SOPHYA::abs(expsignal);
double mine, maxe;
......@@ -691,15 +696,15 @@ int ACxSetFitter::doCxfit(string outfilenamecx, string outcheckfilenamecx, bool
ofr <<setw(4)<<ii+1<<" ERROR FIT RC="<<rcfit<<" nstep="<<mFit.GetNStep()<<endl;
if (_prtlevel_>0) mFit.PrintFitErr(rcfit);
}
if (pox) {
/* DEL if (pox) {
for(size_t j=0; j<NTRK; j++) {
TVector< complex<double> > expsignal = gxi2.getExpectedSignal(j, v_phase[ii], v_Acx[ii][j]);
sprintf(oname,"simcx_%d_%d",(int)ii+1,(int)j+1);
(*pox) << PPFNameTag(oname)<<expsignal;
}
}
} */
}
if (pox) delete pox;
// DEL if (pox) delete pox;
fit_cx_done=true;
acxd_.v_cxbeams=v_cxbeams;
acxd_.v_phase=v_phase;
......@@ -708,10 +713,40 @@ int ACxSetFitter::doCxfit(string outfilenamecx, string outcheckfilenamecx, bool
return 0;
}
int ACxSetFitter::saveExpectedCx(string outcheckfilename)
{
cout << "ACxSetFitter::saveExpectedCx() saving expected cross-cor (and visi-data) to file "<<outcheckfilename<<endl;
POutPersist pox(outcheckfilename);
size_t NB_CXCORS=acxd_.getNbCrossCor();
size_t NTRK = acxd_.NbTrk();
char oname[32];
for(size_t ii=0; ii<NB_CXCORS; ii++) {
CxBeam cxbeam=v_cxbeams[ii];
MyCxSignal cxsig( acxd_.v_time_data, acxd_.vv_cxdata, acxd_.vv_cxerr,
tks_.v_interp_elev, tks_.v_interp_sazim, cxbeam, ii);
for(size_t j=0; j<NTRK; j++) {
TVector< complex<double> > signal = cxsig.getDataSignal(j);
sprintf(oname,"cx_%d_%d",(int)ii+1,(int)j+1);
pox << PPFNameTag(oname)<<signal;
TVector< complex<double> > expsignal = cxsig.getExpectedSignal(j, v_phase[ii], v_Acx[ii][j]);
sprintf(oname,"simcx_%d_%d",(int)ii+1,(int)j+1);
pox << PPFNameTag(oname)<<expsignal;
if (ii==0) {
Vector tmvec = cxsig.getTimeVec(j);
sprintf(oname,"tim_%d",(int)j+1);
pox << PPFNameTag(oname)<<tmvec;
}
}
}
return 0;
}
//------------------------ CxBaselineFitter -------------------------------------
CxBaselineFitter::CxBaselineFitter(vector<AcxDataSet> & v_data, vector<TrackSet> & v_tks)
: v_acxd(v_data), v_trks(v_tks), tot_ntrks(0), fit_done(false), xi2red(-9.e9)
: v_acxd(v_data), v_trks(v_tks), tot_ntrks(0), fit_done(false), xi2red(-9.e9),
bestfitparam(NULL), err_bestfitparam(NULL)
{
if (v_acxd.size() != v_trks.size())
throw ParmError("CxBaselineFitter::CxBaselineFitter(v_data, v_tks) NOT same size v_data,v_tks ");
......@@ -724,21 +759,24 @@ CxBaselineFitter::CxBaselineFitter(vector<AcxDataSet> & v_data, vector<TrackSet>
if (tot_ntrks<1)
throw ParmError("CxBaselineFitter::CxBaselineFitter(v_data, v_tks) 0 tracks ! tot_ntrks<1 ");
size_t nparam = 4*(v_acxd[0].getNbAutoCor()-1);
bestfitparam = new double[nparam];
err_bestfitparam = new double[nparam];
}
CxBaselineFitter::~CxBaselineFitter()
{
if (bestfitparam) delete[] bestfitparam;
if (err_bestfitparam) delete[] err_bestfitparam;
}
int CxBaselineFitter::dofit(string outfilename, string outcheckfilename)
int CxBaselineFitter::dofit(string outfilename)
{
size_t NB_ANTENNES=v_acxd[0].getNbAutoCor(); // nombre d'antennes
size_t NB_CXCORS=v_acxd[0].getNbCrossCor();
cout << "======================================================================================"<<endl;
cout << "------- CxBaselineFitter::dofit() Performing baseline/phase fit on the 6 cross-cors "<<" TotNbTracks="<<tot_ntrks<<endl;
POutPersist * pox = NULL ;
if (outcheckfilename.length()>0) {
cout << "... expected cross-cor for fitted params will be saved to file "<<outcheckfilename<<endl;
pox = new POutPersist(outcheckfilename);
}
ofstream ofr(outfilename.c_str());
ofr << "#### Fitted phases and baseline-shifts (CxBaselineFitter::dofit() ) "<<endl
<< "## NumAntenna Phase BaselineShiftX BaselineShiftY BaselineShiftZ (Phase in degree, BaselineShift in meter) "<<endl;
......@@ -772,11 +810,17 @@ int CxBaselineFitter::dofit(string outfilename, string outcheckfilename)
if (NB_ANTENNES != 4)
throw PError("CxBaselineFitter::dofit() NB_ANTENNES != 4 Current version works only for 4 antenna");
v_phases.resize(NB_ANTENNES-1);
v_err_phases.resize(NB_ANTENNES-1);
v_baselineshits.resize(NB_ANTENNES-1);
v_err_baselineshits.resize(NB_ANTENNES-1);
// SetParam(int n,double value, double step,double min=1., double max=-1.);
for(size_t i=0; i<(NB_ANTENNES-1); i++) {
char pname[32];
sprintf(pname,"Phase_%d",(int)(i+2));
mFit.SetParam(i,pname,v_acxd[0].v_phase[i],M_PI/180.,0.,2.5*M_PI);
v_phases[i]=v_acxd[0].v_phase[i]; v_err_phases[i]=0.;
sprintf(pname,"BaselineShift_X_%d",(int)(i+2));
mFit.SetParam(3+3*i,pname,0.,0.02,-0.25,0.25);
sprintf(pname,"BaselineShift_Y_%d",(int)(i+2));
......@@ -790,7 +834,52 @@ int CxBaselineFitter::dofit(string outfilename, string outcheckfilename)
cout<< "------ Fit result Reduce_Chisquare = " << mFit.GetChi2Red()<< " nstep="<<mFit.GetNStep() << " rc="<<rcfit<<endl;
mFit.PrintFit();
if (pox) delete pox;
for(size_t j=0; j<4; j++)
for(size_t i=0; i<(NB_ANTENNES-1); i++) {
bestfitparam[j*3+i]=mFit.GetParmErr(j*3+i);
err_bestfitparam[j*3+i]=mFit.GetParmErr(j*3+i);
}
for(size_t i=0; i<(NB_ANTENNES-1); i++) {
v_phases[i]=mFit.GetParm(i);
v_err_phases[i]=mFit.GetParmErr(i);
double xs=mFit.GetParm(i*3+3);
double exs=mFit.GetParmErr(i*3+3);
double ys=mFit.GetParm(i*3+4);
double eys=mFit.GetParmErr(i*3+4);
double zs=mFit.GetParm(i*3+5);
double ezs=mFit.GetParmErr(i*3+5);
v_baselineshits[i]=Vector3d(xs,ys,zs);
v_err_baselineshits[i]=Vector3d(exs,eys,ezs);
}
fit_done=true;
return 0;
}
int CxBaselineFitter::saveExpectedCx(string outcheckfilename)
{
cout << "CxBaselineFitter::saveExpectedCx() saving expected cross-cor (and visi-data) to file "<<outcheckfilename<<endl;
POutPersist pox(outcheckfilename);
size_t NB_CXCORS=v_acxd[0].getNbCrossCor();
char oname[48];
My6CxSignalsB cxsigb(v_acxd, v_trks);
for(size_t i=0; i<v_acxd.size(); i++) {
for(size_t j=0; j<v_acxd[i].NbTrk(); j++) {
Vector tmvec = cxsigb.getTimeVec(i,j);
sprintf(oname,"tim_%d_%d",(int)i+1,(int)j+1);
pox << PPFNameTag(oname)<<tmvec;
for(size_t kcx=0; kcx<NB_CXCORS; kcx++) {
TVector< complex<double> > signal = cxsigb.getDataSignal(i,j,kcx);
sprintf(oname,"cx_%d_%d_%d",(int)kcx+1,(int)i+1,(int)j+1);
pox << PPFNameTag(oname)<<signal;
TVector< complex<double> > expsignal = cxsigb.getExpectedSignal(i,j,kcx,bestfitparam);
sprintf(oname,"simcx_%d_%d_%d",(int)kcx+1,(int)i+1,(int)j+1);
pox << PPFNameTag(oname)<<expsignal;
}
}
}
return 0;
}
......@@ -128,10 +128,11 @@ class ACxSetFitter {
public:
ACxSetFitter(AcxDataSet & data, TrackSet & tks);
int doACfit(string outfilename); // outfilename : fit results
int doCxfit(string outfilenamecx, string outcheckfilenamecx, bool useAac=true); // outfilename : fit results
int doCxfit(string outfilenamecx, bool useAac=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
int saveExpectedAC(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; }
......@@ -170,8 +171,11 @@ public:
class CxBaselineFitter {
public:
CxBaselineFitter(vector<AcxDataSet> & v_data, vector<TrackSet> & v_tks);
~CxBaselineFitter();
int dofit(string outfilename);
int dofit(string outfilename, string outcheckfilename);
int saveExpectedCx(string outcheckfilename); // outfilename PPF (or fits ?) file with the expected signals
vector<AcxDataSet> & v_acxd;
vector<TrackSet> & v_trks;
......@@ -181,7 +185,11 @@ public:
int rcfit;
double xi2red;
vector<double> v_phases;
vector<double> v_err_phases;
vector< Vector3d > v_baselineshits;
vector< Vector3d > v_err_baselineshits;
double * bestfitparam;
double * err_bestfitparam;
};
//-----------------------------------------------------------------------
......
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