Skip to content

GitLab

Commit 1bbe957d authored by Reza ANSARI's avatar Reza ANSARI
Browse files

Implementation option fix-baseline ou fix-baseline-XY pour le fit combine des... 

Implementation option fix-baseline ou fix-baseline-XY pour le fit combine des 6 cross-corr et fitSimplex, Reza 05/03/2019
parent 7c4759cc
Hide whitespace changes
Inline Side-by-side
Showing with 87 additions and 55 deletions
gcxfitbaseline.h
View file @ 1bbe957d
......@@ -16,8 +16,8 @@
class My6CxSignalsB {
public:
My6CxSignalsB(vector< AcxDataSet> & acxd, vector< TrackSet > & trks, bool fgno_aphi=false)
: acxd_(acxd), trk_(trks), prtlevel_(0), phlinfac_(1./250.), fgno_aphi_(fgno_aphi)
My6CxSignalsB(vector< AcxDataSet> & acxd, vector< TrackSet > & trks, bool fg_fix_xy=false)
: acxd_(acxd), trk_(trks), prtlevel_(0), phlinfac_(1./250.), fg_fix_xy_(fg_fix_xy)
{
if (acxd_.size() != trk_.size())
throw ParmError("My6CxSignalsB::My6CxSignalsB(acxd,trks)/ERROR acxd_.size() != trk_.size()");
......@@ -146,29 +146,18 @@ public:
size_t nb0[6]={0,1,2,1,2,2};
size_t ja=0, jb=0;
if (fgno_aphi_) {
if (kcx<3) {
phi0=parm[kcx]; aphi=0.;
baseline_shift=Vector3d(parm[3*kcx+3], parm[3*kcx+4], parm[3*kcx+5]);
}
else {
ja=na0[kcx] , jb=nb0[kcx];
phi0=parm[jb]-parm[ja]; aphi=0.;
baseline_shift=Vector3d(parm[3*jb+3]-parm[3*ja+3], parm[3*jb+4]-parm[3*ja+4], parm[3*jb+5]-parm[3*ja+5]);
}
if (kcx<3) {
phi0=parm[2*kcx];
aphi=parm[2*kcx+1];
if (fg_fix_xy_) baseline_shift=Vector3d(0., 0., parm[3*kcx+6]);
else baseline_shift=Vector3d(parm[3*kcx+6], parm[3*kcx+7], parm[3*kcx+8]);
}
else {
if (kcx<3) {
phi0=parm[2*kcx];
aphi=parm[2*kcx+1];
baseline_shift=Vector3d(parm[3*kcx+6], parm[3*kcx+7], parm[3*kcx+8]);
}
else {
ja=na0[kcx] , jb=nb0[kcx];
phi0=parm[2*jb]-parm[2*ja];
aphi=parm[2*jb+1]-parm[2*ja+1];
baseline_shift=Vector3d(parm[3*jb+6]-parm[3*ja+6], parm[3*jb+7]-parm[3*ja+7], parm[3*jb+8]-parm[3*ja+8]);
}
ja=na0[kcx]; jb=nb0[kcx];
phi0=parm[2*jb]-parm[2*ja];
aphi=parm[2*jb+1]-parm[2*ja+1];
if (fg_fix_xy_) baseline_shift=Vector3d(0., 0., parm[3*jb+6]-parm[3*ja+6]);
else baseline_shift=Vector3d(parm[3*jb+6]-parm[3*ja+6], parm[3*jb+7]-parm[3*ja+7], parm[3*jb+8]-parm[3*ja+8]);
}
if (prtlevel_>1) {
cout << "My6CxSignalsB::getFromFitParam(kcx="<<kcx<<") ja="<<ja<<" jb="<<jb
......@@ -254,7 +243,7 @@ public:
vector< size_t > v_npts_xi2;
int prtlevel_;
double phlinfac_;
bool fgno_aphi_; // if true -> pas de terme a_phi (variation lineaire de phase avec al frequence ds le tableau param de getXi2
bool fg_fix_xy_; // if true -> X , Y dans les baselines (seul z varie)
size_t ndatapoints_;
};
......@@ -298,9 +287,9 @@ public:
//-----------------------------------------------------------------
class My6CxMinZFunc : public MinZFunction, public My6CxSignalsB {
public:
My6CxMinZFunc(vector< AcxDataSet> & acxd, vector< TrackSet > & trks, bool fgno_aphi=false)
: MinZFunction((fgno_aphi?12:15)) , // 12 = 3 (phi0) + 3*3 ; 15 = 3 * 2 (phi0,aphi) + 3*3 position shifts
My6CxSignalsB(acxd, trks, fgno_aphi)
My6CxMinZFunc(vector< AcxDataSet> & acxd, vector< TrackSet > & trks, bool fg_fix_xy=false)
: MinZFunction((fg_fix_xy?9:15)) , // 9 = 3*(2 + 1) ( phi0,aphi,z ) ; 15 = 3 * 2 (phi0,aphi) + 3*3 position shifts
My6CxSignalsB(acxd, trks, fg_fix_xy)
{
}
virtual double Value(double const parm[])
......
trkacxfit.cc
View file @ 1bbe957d
......@@ -96,7 +96,8 @@ static bool fg_phi0only = true; // if false, linear varying phase with frequ
// ---- fit simultane sur les 6 cross-correlations
static bool do_baselinefit = false; // if true , perform baseline fit
static bool fg_fixbaseline = false; // if true , do phases fit with fixed baselines
static bool fg_fixbaseline = false; // if true , do phases fit with fixed baselines
static bool fg_fix_xy = false; // if true , do phases fit with z-baselines free only (fixed x,y-baselines)
static bool do_baselinesimplex = false; // if true , perform baseline determination using Simplex minimisation
//--- End of list of global parameters
......@@ -151,17 +152,20 @@ int main (int narg, char* arg[])
string cxbofile = "cxb6_"+outfilename;
string cxbckfile = "cxb6_"+checkfilename;
CxBaselineFitter cxbfit(v_acxd, v_trk);
int fgfixb=0;
if (fg_fixbaseline) fgfixb=2;
else if (fg_fix_xy) fgfixb=1;
//cxbfit.doSimplexMinimize();
cxbfit.initFitParams();
cxbfit.saveExpectedCx(cxbckfile);
if (do_baselinesimplex) {
cxbfit.doSimplexMinimize();
cxbfit.doSimplexMinimize(fgfixb);
cxbckfile="cxb6s_"+checkfilename;
cxbfit.saveExpectedCx(cxbckfile);
// cxbfit.doCheck();
}
if (do_baselinefit) {
cxbfit.dofit(cxbofile,fg_fixbaseline,fg_phi0only);
cxbfit.dofit(cxbofile,fgfixb,fg_phi0only);
cxbckfile = "cxb6f_"+checkfilename;
cxbfit.saveExpectedCx(cxbckfile);
}
......@@ -236,6 +240,7 @@ int decode_args(int narg, char* arg[])
<<" Linear phi model : Phi(freq)=Phi_0 + a_Phi*(freq-1250)/250 \n"
<< " -docx6f : Perform simultaneous fit over 6 cross-corr (baseline and phases determination) - default NO \n"
<< " -fixb : perform the previous fit (docx6f) with fixed baselines - default NO \n"
<< " -fixbxy : perform the previous fit (docx6f) with fixed z-baseline free only (fixed x,y) - default NO \n"
<< " -docx6s : try simplex minimisation over 6 cross-corr (baseline and phases determination) - default NO \n"
<< " -ngb : Use non gaussian beam profile (Bessel j1(angle) - default Gaussian beam) \n"
<< " -D dish_diameter : define effective dish diameter (in meter D*eff , def=4.5) \n"
......@@ -254,6 +259,7 @@ int decode_args(int narg, char* arg[])
fggaussbeam=true;
do_baselinefit=false;
fg_fixbaseline=false;
fg_fix_xy = false;
do_baselinesimplex=false;
fg_cxB0=true;
......@@ -282,6 +288,9 @@ int decode_args(int narg, char* arg[])
else if (fbo=="-fixb") { // Perform the previous docx6f fit with fixed baselines (phases only)
fg_fixbaseline=true; arg++; narg--; lastargs.clear();
}
else if (fbo=="-fixbxy") { // Perform the previous docx6f fit with fixed baselines (phases only)
fg_fix_xy=true; arg++; narg--; lastargs.clear();
}
else if (fbo=="-docx6s") { // Perform simultaneous fit over 6 cross-corr
do_baselinesimplex=true; arg++; narg--; lastargs.clear();
}
......@@ -318,8 +327,10 @@ int decode_args(int narg, char* arg[])
<< " Phase model : Phase(freq)= " << (fg_phi0only?"Phi0":"Phi0+aPhi*(freq-1250)/250")
<< " Fit B0 ? "<<(fg_cxB0?"No (B=(0,0)":"Yes")<<endl;
cout << " Perform baseline fit ? " << (do_baselinefit?"Yes":"No")
<<(fg_fixbaseline?" With FIXED baselines":" (Phases & baselines)")
<<" Simplex-Minimisation ? "<<(do_baselinesimplex?"Yes":"No")<<endl;
<<" Simplex-Minimisation ? "<<(do_baselinesimplex?"Yes":"No");
if (fg_fixbaseline) cout << " With FIXED baselines"<<endl;
else if (fg_fix_xy) cout << " With FIXED X,Y baselines"<<endl;
else cout << " fit Phases & baselines "<<endl;
cout << " --- TrackSetFiles: (NbFiles="<<lastargs.size()<<" default directory: "<<default_input_path<<")"<<endl;
trksetfiles = lastargs;
for (size_t i=0; i<trksetfiles.size(); i++) {
......
trkfit.cc
View file @ 1bbe957d
......@@ -957,6 +957,8 @@ CxBaselineFitter::~CxBaselineFitter()
void CxBaselineFitter::initFitParams()
{
// Positions en z des 3 antennes 2,3,4 / 1 , en metres
double z_antennes[3]={-0.10,+0.05,-0.02};
//DBG cout << " *DBG* CxBaselineFitter::initFitParams() v_acxd[0].v_phase.size()="<<v_acxd[0].v_phase.size()<<endl;
size_t NB_ANTENNES=v_acxd[0].getNbAutoCor(); // nombre d'antennes
size_t NB_CXCORS=v_acxd[0].getNbCrossCor();
......@@ -985,7 +987,7 @@ void CxBaselineFitter::initFitParams()
}
int CxBaselineFitter::dofit(string outfilename, bool fgfixbaseline, bool fgphi0only)
int CxBaselineFitter::dofit(string outfilename, int fgfixbaseline, bool fgphi0only)
{
size_t NB_ANTENNES=v_acxd[0].getNbAutoCor(); // nombre d'antennes
size_t NB_CXCORS=v_acxd[0].getNbCrossCor();
......@@ -1038,11 +1040,18 @@ int CxBaselineFitter::dofit(string outfilename, bool fgfixbaseline, bool fgphi0o
mFit.SetParam(7+3*i,pname,v_baselineshits[i].Y(),0.01,-0.25,0.25);
sprintf(pname,"BaselineShift_Z_%d",(int)(i+2));
mFit.SetParam(8+3*i,pname,v_baselineshits[i].Z(),0.01,-0.25,0.25);
if (fgfixbaseline) {
if (fgfixbaseline==1) {
cout << " ... fitting Z-baseline and phases , fixed X,Y baselines "<<endl;
mFit.SetFix(6+3*i); mFit.SetFix(7+3*i);
}
else if (fgfixbaseline==2) {
cout << " ... fitting phases only, fixed baselines "<<endl;
mFit.SetFix(6+3*i); mFit.SetFix(7+3*i); mFit.SetFix(8+3*i);
}
}
int ndataused=0;
double xi2redini=gxi2.getXi2(bestfitparam, ndataused);
cout << " Xi2Red value for the initial parameters guess= " <<xi2redini/(double)ndataused<<endl;
cout << " Performing the fit (tot_npoints_fit= "<<tot_npoints_fit<<" ?= (npoints2="<<npoints2<<") ..."<< endl;
rcfit = mFit.doFit(); xi2red=-99999.;
cout<< "------ Fit result Reduce_Chisquare = " << mFit.GetChi2Red()<< " nstep="<<mFit.GetNStep() << " rc="<<rcfit<<endl;
......@@ -1081,7 +1090,7 @@ int CxBaselineFitter::dofit(string outfilename, bool fgfixbaseline, bool fgphi0o
return 0;
}
int CxBaselineFitter::doSimplexMinimize()
int CxBaselineFitter::doSimplexMinimize(int fg_fixebaseline)
{
size_t NB_ANTENNES=v_acxd[0].getNbAutoCor(); // nombre d'antennes
size_t NB_CXCORS=v_acxd[0].getNbCrossCor();
......@@ -1091,18 +1100,26 @@ int CxBaselineFitter::doSimplexMinimize()
if (NB_ANTENNES != 4)
throw PError("CxBaselineFitter::doSimplexMinimize() NB_ANTENNES != 4 Current version works only for 4 antenna");
My6CxMinZFunc mzfunc(v_acxd, v_trks, true);
bool fgfixbasel=false;
if (fg_fixebaseline>0) fgfixbasel=true;
My6CxMinZFunc mzfunc(v_acxd, v_trks, fgfixbasel);
MinZSimplex simplex(&mzfunc);
// Guess the center and step for constructing the initial simplex
size_t nparam = 4*(NB_ANTENNES-1);
size_t nparam = 5*(NB_ANTENNES-1);
size_t npar_pas = 5;
if (fgfixbasel) { nparam = 3*(NB_ANTENNES-1); npar_pas=3; }
Vector P0(nparam);
Vector step(nparam);
for(size_t i=0; i<(NB_ANTENNES-1); i++) {
P0(i)=v_acxd[0].v_phase[i];
step(i)=M_PI/6.;
for(size_t j=0;j<3;j++) {
P0((i+1)*3+j)=0.;
step((i+1)*3+j)=0.05;
P0(2*i)=v_phi_0[i];
step(2*i)=M_PI/20.;
P0(2*i+1)=v_a_phi[i];
step(2*i+1)=M_PI/30.;
if (fgfixbasel) {
for(size_t j=0; j<3; j++) { P0(i+6+j)=0.; step(i+6+j)=0.02; }
}
else {
for(size_t j=0; j<9; j++) { P0(i+6+j)=0.; step(i+6+j)=0.02; }
}
}
cout << " Initial Point: "<<P0.Transpose()<<endl;
......@@ -1123,19 +1140,28 @@ int CxBaselineFitter::doSimplexMinimize()
cout << " Converged: NStep= " << simplex.NbIter() << " Best Xi2="<< mzfunc.Value(oparm.Data()) << endl;
simplex_done=true;
for(size_t i=0; i<(NB_ANTENNES-1); i++) {
v_phi_0[i]=oparm(i);
double xs=oparm(i*3+3);
double ys=oparm(i*3+4);
double zs=oparm(i*3+5);
v_phi_0[i]=oparm(2*i);
v_a_phi[i]=oparm(2*i+1);
double xs=0., ys=0., zs=0.;
if (fgfixbasel) {
xs=0.; ys=0.;
zs=oparm(6+i);
}
else {
xs=oparm(6+3*i);
ys=oparm(7+3*i);
zs=oparm(8+i);
}
v_baselineshits[i]=Vector3d(xs,ys,zs);
cout << " ANTENNE["<<i+2<<"] : Phase="<<v_phi_0[i]<<" BaseLineShift="<<v_baselineshits[i]<<endl;
cout << " ANTENNE["<<i+2<<"] : Phi0="<<v_phi_0[i]<<" aPhi="<<v_a_phi[i]<<" BaseLineShift="<<v_baselineshits[i]<<endl;
}
}
return 0;
}
int CxBaselineFitter::doCheck()
/* A SUPPRIMER
int CxBaselineFitter::doCheck(int fg_fixebaseline)
{
size_t NB_ANTENNES=v_acxd[0].getNbAutoCor(); // nombre d'antennes
size_t NB_CXCORS=v_acxd[0].getNbCrossCor();
......@@ -1145,14 +1171,18 @@ int CxBaselineFitter::doCheck()
if (NB_ANTENNES != 4)
throw PError("CxBaselineFitter::doCheck() NB_ANTENNES != 4 Current version works only for 4 antenna");
My6CxMinZFunc mzfunc(v_acxd, v_trks, true); // true : Pas de aphi ds les tableaux param
bool fgfixbasel=false;
if (fg_fixebaseline>0) fgfixbasel=true;
My6CxMinZFunc mzfunc(v_acxd, v_trks, fgfixbasel); // true : Pas de aphi ds les tableaux param
mzfunc.SetPrintLevel(_prtlevel_);
// Guess the center and step for constructing the initial simplex
size_t nparam = 4*(NB_ANTENNES-1);
size_t nparam = 5*(NB_ANTENNES-1);
if (fgfixbasel) nparam = 3*(NB_ANTENNES-1);
Vector P0(nparam), PC(nparam);
Vector step(nparam);
for(size_t i=0; i<(NB_ANTENNES-1); i++) {
P0(i)=v_acxd[0].v_phase[i];
P0(i)=v_phi_0[i];
P0(1)=v_a_phi[i];
step(i)=M_PI;
for(size_t j=0;j<3;j++) {
P0((i+1)*3+j)=0.;
......@@ -1197,6 +1227,7 @@ int CxBaselineFitter::doCheck()
return 0;
}
*/
int CxBaselineFitter::saveExpectedCx(string outcheckfilename)
{
......
trkfit.h
View file @ 1bbe957d
......@@ -205,11 +205,12 @@ public:
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);
// if fgfixbaseline =0 -> free baselines (fit x,y,z) = 1 -> Fix X,Y-baselines =2 -> fix baselines (phases only fit)
// if fgphi0only : phases independent of frequency , aphi=0
int dofit(string outfilename, int fgfixebaseline=0, bool fgphi0only=true);
int doCheck();
int doSimplexMinimize();
int doCheck(int fg_fixebaseline=0);
int doSimplexMinimize(int fg_fixebaseline=0);
int saveExpectedCx(string outcheckfilename); // outfilename PPF (or fits ?) file with the expected signals
......
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