Skip to content

GitLab

Commit 2511dec1 authored by Reza ANSARI's avatar Reza ANSARI
Browse files

Suite implementation fits multi-frequence et decodage d'un argument optionnel... 

Suite implementation fits multi-frequence et decodage d'un argument optionnel NOAC, NOCX NOACCX du datacard @trk, Reza 25/02/2019
parent 36678b83
Hide whitespace changes
Inline Side-by-side
Showing with 215 additions and 99 deletions
gacfit.h
View file @ 2511dec1
......@@ -17,10 +17,10 @@ static int debug_gacfit=0;
class MyACSignal {
public:
MyACSignal(vector< vector<double> > & v_time_data_, vector< vector< vector<double> > > & vv_data_,
vector< vector< vector<double> > > & vv_err_, vector<double> & v_freqs_,
vector< vector< vector<double> > > & vv_err_, vector<double> & v_freqs_, vector<bool> & v_noAC_,
vector< SLinInterp1D > & v_interp_elev_, vector< SLinInterp1D > & v_interp_sazim_,
size_t nac=0, bool fggauss_beam=true)
: v_time_data(v_time_data_), vv_data(vv_data_), vv_err(vv_err_), v_freqs(v_freqs_),
: v_time_data(v_time_data_), vv_data(vv_data_), vv_err(vv_err_), v_freqs(v_freqs_), v_noAC(v_noAC_),
v_interp_elev(v_interp_elev_), v_interp_sazim(v_interp_sazim_),
nac_(nac), fggauss_beam_(fggauss_beam)
{
......@@ -117,6 +117,7 @@ public:
double xi2=0.;
npts=0;
for(size_t j=0; j<v_time_data.size(); j++) {
if (v_noAC[j]) continue;
vector<double> sig;
npts += getExpectedSignal(j, sig, parm);
vector<double> & vdata = vv_data[j][nac_];
......@@ -133,7 +134,8 @@ public:
vector< vector< vector<double> > > & vv_data;
vector< vector< vector<double> > > & vv_err;
vector<double> & v_freqs;
vector<double> v_lambdas;
vector<double> v_lambdas;
vector<bool> & v_noAC;
vector< SLinInterp1D > & v_interp_elev;
vector< SLinInterp1D > & v_interp_sazim;
size_t nac_;
......@@ -144,11 +146,12 @@ public:
class MyACGenXi2 : public GeneralXi2, MyACSignal {
public:
MyACGenXi2(vector< vector<double> > & v_time_data_, vector< vector< vector<double> > > & vv_data_,
vector< vector< vector<double> > > & vv_err_, vector<double> & v_freqs_,
vector< vector< vector<double> > > & vv_err_, vector<double> & v_freqs_, vector<bool> & v_noAC_,
vector< SLinInterp1D > & v_interp_elev_, vector< SLinInterp1D > & v_interp_sazim_,
size_t nac=0, bool fggauss_beam=true)
: GeneralXi2(3+2*v_time_data_.size()) ,
MyACSignal(v_time_data_, vv_data_, vv_err_, v_freqs_, v_interp_elev_, v_interp_sazim_, nac, fggauss_beam)
MyACSignal(v_time_data_, vv_data_, vv_err_, v_freqs_, v_noAC_,
v_interp_elev_, v_interp_sazim_, nac, fggauss_beam)
{
}
......@@ -169,11 +172,12 @@ public:
class MyACMinZ : public MinZFunction, public MyACSignal {
public:
MyACMinZ(vector< vector<double> > & v_time_data_, vector< vector< vector<double> > > & vv_data_,
vector< vector< vector<double> > > & vv_err_, vector<double> & v_freqs_,
vector< SLinInterp1D > & v_interp_elev_, vector< SLinInterp1D > & v_interp_sazim_,
size_t nac=0, bool fggauss_beam=true)
vector< vector< vector<double> > > & vv_err_, vector<double> & v_freqs_, vector<bool> & v_noAC_,
vector< SLinInterp1D > & v_interp_elev_, vector< SLinInterp1D > & v_interp_sazim_,
size_t nac=0, bool fggauss_beam=true)
: MinZFunction(3+2*v_time_data_.size()) ,
MyACSignal(v_time_data_, vv_data_, vv_err_, v_freqs_, v_interp_elev_, v_interp_sazim_, nac, fggauss_beam)
MyACSignal(v_time_data_, vv_data_, vv_err_, v_freqs_, v_noAC_,
v_interp_elev_, v_interp_sazim_, nac, fggauss_beam)
{
}
virtual double Value(double const parm[])
......@@ -184,4 +188,44 @@ public:
};
#ifdef TKF_AVEC_MINUIT
// Ajustement avec Minuit
class MyAC_Xi2_Minuit : public ROOT::Minuit2::FCNBase , public MyACSignal {
public:
MyAC_Xi2_Minuit(vector< vector<double> > & v_time_data_, vector< vector< vector<double> > > & vv_data_,
vector< vector< vector<double> > > & vv_err_, vector<double> & v_freqs_, vector<bool> & v_noAC_,
vector< SLinInterp1D > & v_interp_elev_, vector< SLinInterp1D > & v_interp_sazim_,
size_t nac=0, bool fggauss_beam=true)
: MyACSignal(v_time_data_, vv_data_, vv_err_, v_freqs_, v_noAC_,
v_interp_elev_, v_interp_sazim_, nac, fggauss_beam)
{
nbparam_=3+2*v_time_data_.size();
myparam_ = new double[nbparam_];
}
virtual ~MyAC_Xi2_Minuit()
{
delete[] myparam_;
}
/*
param[0] = D (Dish diameter)
parm[1] = delta_elev = theta
parm[2] = azim = phi
parm[3+2*j] = Aj (Amplitude) for track j
parm[4+2*j] = Bj (Baseline / Offset) for track j
*/
virtual double operator() (const std::vector< double > & p) const
{
if (p.size() != nbparam_) throw ParmError("MyAC_Xi2_Minuit::operator()(p) p.size() != nbparam_");
for(size_t i=0; i<nbparam_; i++) myparam_[i]=p[i];
int ndataused=0;
return getXi2(myparam_, ndataused);
}
size_t nbparam_;
double * myparam_;
};
#endif /* TKF_AVEC_MINUIT */
#endif
gcxfit.h
View file @ 2511dec1
......@@ -49,7 +49,8 @@ public:
return rvs;
}
// Linearly varying phase with frequency
//---- This method should be identical to the one in My6CxSignalsB (gcxfitbaseline.h)
// Model for linearly varying phase with frequency : Phase = Phi0 + a_phi * (freq-1250)/250
inline double getPhase4Freq(double phi0, double a_phi, double freq)
{
return (phi0+a_phi*(freq-1250.)*phlinfac_);
......
gcxfitbaseline.h
View file @ 2511dec1
......@@ -14,8 +14,8 @@
class My6CxSignalsB {
public:
My6CxSignalsB(vector< AcxDataSet> & acxd, vector< TrackSet > & trks)
: acxd_(acxd), trk_(trks), prtlevel_(0)
My6CxSignalsB(vector< AcxDataSet> & acxd, vector< TrackSet > & trks, bool fgno_aphi=false)
: acxd_(acxd), trk_(trks), prtlevel_(0), phlinfac_(1./250.), fgno_aphi_(fgno_aphi)
{
if (acxd_.size() != trk_.size())
throw ParmError("My6CxSignalsB::My6CxSignalsB(acxd,trks)/ERROR acxd_.size() != trk_.size()");
......@@ -68,8 +68,16 @@ public:
return acxd_[i].v_cxbeams[kcx];
}
//---- This method should be identical to the one in MyCxSignal (gcxfit.h)
// Model for linearly varying phase with frequency : Phase = Phi0 + a_phi * (freq-1250)/250
inline double getPhase4Freq(double phi0, double a_phi, double freq)
{
return (phi0+a_phi*(freq-1250.)*phlinfac_);
}
// 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)
virtual int getExpectedSignal(size_t i, size_t j, size_t kcx, vector< complex<double> > & sig,
double phi0, double aphi, Vector3d & baseline_shift)
{
//DBG cout<<"*DBG*A*My6CxSignalsB::getExpectedSignal() j="<<j<<" size="<<v_time_data[j].size()<<" kcx="<<kcx<<endl;
vector<double> & vtime = acxd_[i].v_time_data[j];
......@@ -81,6 +89,9 @@ public:
CxBeam beam = getBeam(i, kcx);
beam.ShiftBaseline(baseline_shift);
beam.setFreq(acxd_[i].v_freqs[j]);
double phase = getPhase4Freq(phi0, aphi, acxd_[i].v_freqs[j]);
if (prtlevel_>1) {
cout << "My6CxSignalsB::getExpectedSignal() A="<<A<<" beam.baseline="<<beam.baseline_
<< " a1.D="<<beam.a1_.getDdish()<<endl;
......@@ -104,21 +115,25 @@ public:
}
/* definition du tableau des parametres de fit
parm[0] = Phase_2
parm[1] = Phase_3
parm[2] = Phase_4
parm[3] = X_shift_2
parm[4] = Y_shift_2
parm[5] = Z_shift_2
parm[6] = X_shift_3
parm[7] = Y_shift_3
parm[8] = Z_shift_3
parm[9] = X_shift_4
parm[10] = Y_shift_4
parm[11] = Z_shift_4
parm[0] = Phi0_2
parm[1] = a_Phi_2
parm[2] = Phi0_3
parm[3] = a_Phi_3
parm[4] = Phi0_4
parm[5] = a_Phi_4
parm[6] = X_shift_2
parm[7] = Y_shift_2
parm[8] = Z_shift_2
parm[9] = X_shift_3
parm[10] = Y_shift_3
parm[11] = Z_shift_3
parm[12] = X_shift_4
parm[13] = Y_shift_4
parm[14] = Z_shift_4
*/
// determine la phase et le baseline_shift a partir du tableau parm, pour la ligne de base kcx (0<=kcx<=5)
void getFromFitParam(size_t kcx, const double* parm, double & phase, Vector3d & baseline_shift)
void getFromFitParam(size_t kcx, const double* parm, double & phi0, double & aphi, Vector3d & baseline_shift)
{
size_t na[6]={0,0,0,1,1,2};
size_t nb[6]={1,2,3,2,3,3};
......@@ -126,19 +141,34 @@ public:
size_t nb0[6]={0,1,2,1,2,2};
size_t ja=0, jb=0;
if (kcx<3) {
phase=parm[kcx];
baseline_shift=Vector3d(parm[3*kcx+3], parm[3*kcx+4], parm[3*kcx+5]);
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]);
}
}
else {
ja=na0[kcx] , jb=nb0[kcx];
phase=parm[jb]-parm[ja];
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];
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]);
}
}
if (prtlevel_>1) {
cout << "My6CxSignalsB::getFromFitParam(kcx="<<kcx<<") ja="<<ja<<" jb="<<jb
<<" ->phase="<<phase<<" BaselineShift="<<baseline_shift<<endl<<" parm[0..11]= ";
for(size_t k=0; k<12; k++) cout << parm[k]<<" ; ";
<<" ->Phi0="<<phi0<<" aPhi="<<aphi<<" BaselineShift="<<baseline_shift<<endl<<" parm[0..11]= ";
for(size_t k=0; k<15; k++) cout << parm[k]<<" ; ";
cout<<endl;
}
......@@ -147,19 +177,19 @@ public:
// 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;
double phi0,aphi;
Vector3d baseline_shift;
getFromFitParam(kcx, parm, phase, baseline_shift);
return getExpectedSignal(i, j, kcx, sig, phase, baseline_shift);
getFromFitParam(kcx, parm, phi0, aphi, baseline_shift);
return getExpectedSignal(i, j, kcx, sig, phi0, aphi, 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)
virtual TVector< complex<double> > getExpectedSignal(size_t i, size_t j, size_t kcx, double phi0, double aphi, Vector3d & baseline_shift)
{
vector< complex<double> > sig;
//DBG cout << "*DBG*getExpectedSignal() : calling getExpectedSignal(j="<<j<<" ,sig,...)"<<endl;
//DBG debug_gacfit=1;
getExpectedSignal(i, j, kcx, sig, phase, baseline_shift);
getExpectedSignal(i, j, kcx, sig, phi0, aphi, baseline_shift);
TVector< complex<double> > rvs(sig);
return rvs;
}
......@@ -180,16 +210,16 @@ public:
double xi2=0.;
npts=0;
double myphase=0.;
double myphi0=0., myaphi=0.;
Vector3d myblshift;
for(size_t kcx=0; kcx<n_cxcor; kcx++) { // loop over the six cross-correlations
v_xi2[kcx]=0.;
double myxi2=0.;
getFromFitParam(kcx, parm, myphase, myblshift);
getFromFitParam(kcx, parm, myphi0, myaphi, myblshift);
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
vector< complex<double> > sig;
npts += getExpectedSignal(i, j, kcx, sig, myphase, myblshift);
npts += getExpectedSignal(i, j, kcx, sig, myphi0, myaphi, myblshift);
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++) {
......@@ -218,30 +248,36 @@ public:
vector< double > v_xi2;
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
};
//-----------------------------------------------------------------
class My6CxGenXi2B : public GeneralXi2, public My6CxSignalsB {
public:
My6CxGenXi2B(vector< AcxDataSet> & acxd, vector< TrackSet > & trks)
: GeneralXi2(12) , // 12 = 3 phases + 3*3 position shifts
: GeneralXi2(15) , // nb_param= 15 = 3 * 2 (phi0,aphi) + 3*3 position shifts
My6CxSignalsB(acxd, trks)
{
}
/* definition du tableau des parametres de fit
parm[0] = Phase_2 Phase Antenne 2 = DeltaPhi_12
parm[1] = Phase_3 Phase Antenne 3 = DeltaPhi_13
parm[2] = Phase_4 Phase Antenne 4 = DeltaPhi_14
parm[3] = X_shift_2 Shift (/nominal) X Antenne 2 ( DeltaX_12 )
parm[4] = Y_shift_2 Shift (/nominal) Y Antenne 2 ( DeltaY_12 )
parm[5] = Z_shift_2 Shift (/nominal) Z Antenne 2 ( DeltaZ_12 )
parm[6] = X_shift_3 Shift (/nominal) X Antenne 3 ( DeltaX_13 )
parm[7] = Y_shift_3 Shift (/nominal) Y Antenne 3 ( DeltaY_13 )
parm[8] = Z_shift_3 Shift (/nominal) Z Antenne 2 ( DeltaZ_13 )
parm[9] = X_shift_4 Shift (/nominal) X Antenne 4 ( DeltaX_14 )
parm[10] = Y_shift_4 Shift (/nominal) Y Antenne 4 ( DeltaY_14 )
parm[11] = Z_shift_4 Shift (/nominal) Z Antenne 4 ( DeltaZ_14 )
parm[0] = Phi0_2 Phi0 Antenne 2
parm[1] = a_Phi_2 aPhi (slope) Antenne 2
parm[2] = Phi0_3 Phi0 Antenne 3
parm[3] = a_Phi_3 aPhi (slope) Antenne 3
parm[4] = Phi0_4 Phi0 Antenne 4
parm[5] = a_Phi_4 aPhi (slope) Antenne 4
parm[6] = X_shift_2 Shift (/nominal) X Antenne 2 ( DeltaX_12 )
parm[7] = Y_shift_2 Shift (/nominal) Y Antenne 2 ( DeltaY_12 )
parm[8] = Z_shift_2 Shift (/nominal) Z Antenne 2 ( DeltaZ_12 )
parm[9] = X_shift_3 Shift (/nominal) X Antenne 3 ( DeltaX_13 )
parm[10] = Y_shift_3 Shift (/nominal) Y Antenne 3 ( DeltaY_13 )
parm[11] = Z_shift_3 Shift (/nominal) Z Antenne 2 ( DeltaZ_13 )
parm[12] = X_shift_4 Shift (/nominal) X Antenne 4 ( DeltaX_14 )
parm[13] = Y_shift_4 Shift (/nominal) Y Antenne 4 ( DeltaY_14 )
parm[14] = Z_shift_4 Shift (/nominal) Z Antenne 4 ( DeltaZ_14 )
*/
virtual double Value(GeneralFitData& data, double* parm, int& ndataused)
......@@ -253,9 +289,9 @@ public:
//-----------------------------------------------------------------
class My6CxMinZFunc : public MinZFunction, public My6CxSignalsB {
public:
My6CxMinZFunc(vector< AcxDataSet> & acxd, vector< TrackSet > & trks)
: MinZFunction(12) , // 12 = 3 phases + 3*3 position shifts
My6CxSignalsB(acxd, trks)
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)
{
}
virtual double Value(double const parm[])
......
trkfit.cc
View file @ 2511dec1
......@@ -43,6 +43,8 @@ static vector<double> * tstart_p_ = NULL;
static vector<double> * tend_p_ = NULL;
static vector<double> * v_freqs_p_ = NULL;
static vector<string> * trkflnm_p_ = NULL;
static vector<bool> * v_noAC_p_ = NULL;
static vector<bool> * v_noCx_p_ = NULL;
static size_t trk_cnt = 0;
static int decode_trkcard(string const& key, string const& toks)
......@@ -56,14 +58,26 @@ static int decode_trkcard(string const& key, string const& toks)
return 1;
}
char flnmdata[256], flnmtrk[256];
char sflags[64];
double ts,te,freq;
sscanf(toks.c_str(),"%s %lg,%lg %lg %s",flnmdata,&ts,&te,&freq,flnmtrk);
sscanf(toks.c_str(),"%s %lg,%lg %lg %s %s",flnmdata,&ts,&te,&freq,flnmtrk,sflags);
dataflnm_p_->push_back(flnmdata);
tstart_p_->push_back(ts*60.);
tend_p_->push_back(te*60.);
v_freqs_p_->push_back(freq);
trkflnm_p_->push_back(flnmtrk);
size_t ll=strlen(sflags);
bool noAC=false;
bool noCx=false;
if (ll>0) {
for(size_t l=0; l<ll; l++) sflags[l]=toupper(sflags[l]);
string sflg=sflags;
if ((sflg == "NOAC")||(sflg=="NOACCX")) noAC=true;
if ((sflg == "NOCX")||(sflg=="NOACCX")) noCx=true;
}
v_noAC_p_->push_back(noAC);
v_noCx_p_->push_back(noCx);
trk_cnt++;
return 0;
}
......@@ -92,9 +106,14 @@ size_t TrkInputDataSet::ReadDatacardFile(string dcfilename)
tend_p_ = &tend;
v_freqs_p_ = &v_freqs;
trkflnm_p_ = &trkflnm;
v_noAC_p_=&v_noAC;
v_noCx_p_=&v_noCx;
trk_cnt = 0;
// @trk visiDataTableFile tstart,tend freq TrackFileName
// @trk visiDataTableFile tstart,tend freq TrackFileName [FLAG]
// tstart , tend in minutes freq in MHz
// optional FLAG = NOAC NOCX NOACCX
// NOAC : don't use for Auto-correlation fit ; NOCX : don't use for cross-cor fits
// NOACCX : don't use for Auto-correlation or cross-cor fits
dc.ReadFile(dcfilename);
if (dc.HasKey("inpath")) { // @inpath InputFilesDirectoryPath
input_base_path = dc.SParam("inpath",0,"");
......@@ -150,7 +169,7 @@ AcxDataSet::AcxDataSet(AcxDataSet const & a)
v_min_data(a.v_min_data), v_max_data(a.v_max_data),
vv_cxdata(a.vv_cxdata), vv_cxerr(a.vv_cxerr),
v_min_cxdata(a.v_min_cxdata), v_max_cxdata(a.v_max_cxdata),
tot_npoints(a.tot_npoints), v_freqs(a.v_freqs),
tot_npoints(a.tot_npoints), v_freqs(a.v_freqs), v_noAC(a.v_noAC), v_noCx(a.v_noCx),
zenang(a.zenang), theta_0(a.theta_0), phi_0(a.phi_0),
v_acbeams(a.v_acbeams), v_cxbeams(a.v_cxbeams),
v_phase(a.v_phase), v_phi_0(a.v_phi_0), v_a_phi(a.v_a_phi), v_Acx(a.v_Acx), v_Bcx(a.v_Bcx)
......@@ -163,7 +182,7 @@ AcxDataSet & AcxDataSet::operator = (AcxDataSet const & a)
v_min_data=a.v_min_data; v_max_data=a.v_max_data;
vv_cxdata=a.vv_cxdata; vv_cxerr=a.vv_cxerr;
v_min_cxdata=a.v_min_cxdata; v_max_cxdata=a.v_max_cxdata;
tot_npoints=a.tot_npoints; v_freqs=a.v_freqs;
tot_npoints=a.tot_npoints; v_freqs=a.v_freqs; v_noAC=a.v_noAC; v_noCx=a.v_noCx;
zenang=a.zenang; theta_0=a.theta_0; phi_0=a.phi_0;
v_acbeams=a.v_acbeams; v_cxbeams=a.v_cxbeams;
v_phase=a.v_phase; v_phi_0=a.v_phi_0; v_a_phi=a.v_a_phi; v_Acx=a.v_Acx; v_Bcx=a.v_Bcx;
......@@ -187,7 +206,7 @@ size_t AcxDataSet::ReadData(TrkInputDataSet & tkds)
v_min_cxdata.resize(tkds.NbTrk());
v_max_cxdata.resize(tkds.NbTrk());
}
v_freqs=tkds.v_freqs;
v_freqs=tkds.v_freqs; v_noAC=tkds.v_noAC; v_noCx=tkds.v_noCx;
zenang=tkds.zenang; theta_0=tkds.theta_0; phi_0=tkds.phi_0;
size_t NB_ANTENNES=getNbAutoCor(); // nombre d'antennes
size_t NB_CXCORS=getNbCrossCor();
......@@ -408,7 +427,7 @@ int ACxSetFitter::doACfit(string outfilename)
gdata.AddData1(acxd_.v_time_data[j][k],v_data[ii][k],v_err[ii][k]); // Fill x, y and error on y
}
}
MyACGenXi2 gxi2( acxd_.v_time_data, acxd_.vv_data, acxd_.vv_err, acxd_.v_freqs,
MyACGenXi2 gxi2( acxd_.v_time_data, acxd_.vv_data, acxd_.vv_err, acxd_.v_freqs, acxd_.v_noAC,
tks_.v_interp_elev, tks_.v_interp_sazim, ii, fggaussbeam_);
GeneralFit mFit(&gxi2);
mFit.SetData(&gdata); // connect data to the fitter , here the data is unused - gxi2 includes its data
......@@ -521,7 +540,7 @@ int ACxSetFitter::saveExpectedAC(string outcheckfilename)
if (_prtlevel_>1)
cout << "... Computing DataSignal & Expected Signal for fitted params and dish "<<ii+1<<endl;
MyACSignal macs(acxd_.v_time_data, acxd_.vv_data, acxd_.vv_err, acxd_.v_freqs,
MyACSignal macs(acxd_.v_time_data, acxd_.vv_data, acxd_.vv_err, acxd_.v_freqs, acxd_.v_noAC,
tks_.v_interp_elev, tks_.v_interp_sazim, ii, fggaussbeam_);
double Ddishfit=v_Ddish[ii];
......@@ -784,14 +803,13 @@ CxBaselineFitter::CxBaselineFitter(vector<AcxDataSet> & v_data, vector<TrackSet>
throw ParmError("CxBaselineFitter::CxBaselineFitter(v_data, v_tks) NOT same size v_data,v_tks ");
if (v_acxd.size() < 1)
throw ParmError("CxBaselineFitter::CxBaselineFitter(v_data, v_tks) v_data.size()<1 ");
v_phases.resize(v_acxd[0].getNbAutoCor()-1);
v_baselineshits.resize(v_acxd[0].getNbAutoCor()-1);
tot_ntrks=0;
for(size_t i=0; i<v_acxd.size(); i++) tot_ntrks+=v_acxd[i].NbTrk();
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);
size_t nparam = 5*(v_acxd[0].getNbAutoCor()-1); // 5 param / antenne , phi0, aphi, dX,dY,dZ
bestfitparam = new double[nparam];
err_bestfitparam = new double[nparam];
......@@ -811,16 +829,21 @@ void CxBaselineFitter::initFitParams()
size_t NB_CXCORS=v_acxd[0].getNbCrossCor();
if (NB_ANTENNES != 4)
throw PError("CxBaselineFitter::initFitParams() NB_ANTENNES != 4 Current version works only for 4 antenna");
v_phases.resize(NB_ANTENNES-1);
v_err_phases.resize(NB_ANTENNES-1);
v_phi_0.resize(v_acxd[0].getNbAutoCor()-1);
v_err_phi_0.resize(NB_ANTENNES-1);
v_a_phi.resize(v_acxd[0].getNbAutoCor()-1);
v_err_a_phi.resize(NB_ANTENNES-1);
v_baselineshits.resize(NB_ANTENNES-1);
v_err_baselineshits.resize(NB_ANTENNES-1);
for(size_t i=0; i<(NB_ANTENNES-1); i++) {
v_phases[i]=v_acxd[0].v_phase[i]; v_err_phases[i]=0.;
v_phi_0[i]=v_acxd[0].v_phi_0[i]; v_err_phi_0[i]=0.;
v_a_phi[i]=v_acxd[0].v_a_phi[i]; v_err_a_phi[i]=0.;
v_baselineshits[i]=Vector3d(0.,0.,0.);
v_err_baselineshits[i]=Vector3d(0.,0.,0.);
bestfitparam[i]=v_phases[i];
err_bestfitparam[i]=0.;
bestfitparam[2*i]=v_phi_0[i];
err_bestfitparam[2*i]=0.;
bestfitparam[2*i+1]=v_a_phi[i];
err_bestfitparam[2*i]=0.;
for(size_t j=0; j<3; j++) {
bestfitparam[3*(i+1)+j]=err_bestfitparam[3*(i+1)+j]=0.;
}
......@@ -829,7 +852,7 @@ void CxBaselineFitter::initFitParams()
}
int CxBaselineFitter::dofit(string outfilename, bool fgfixbaseline)
int CxBaselineFitter::dofit(string outfilename, bool fgfixbaseline, bool fgphi0only)
{
size_t NB_ANTENNES=v_acxd[0].getNbAutoCor(); // nombre d'antennes
size_t NB_CXCORS=v_acxd[0].getNbCrossCor();
......@@ -868,18 +891,21 @@ int CxBaselineFitter::dofit(string outfilename, bool fgfixbaseline)
// 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_phases[i],M_PI/180.,0.,2.5*M_PI);
v_err_phases[i]=0.;
sprintf(pname,"Phi0_%d",(int)(i+2));
mFit.SetParam(2*i,pname,v_phi_0[i],M_PI/180.,0.,2.5*M_PI);
sprintf(pname,"a_Phi_%d",(int)(i+2));
mFit.SetParam(2*i+1,pname,v_a_phi[i],0.1,-15.,15.);
if (fgphi0only) mFit.SetFix(2*i+1, 0.);
v_err_phi_0[i]=0.; v_err_a_phi[i]=0.;
sprintf(pname,"BaselineShift_X_%d",(int)(i+2));
mFit.SetParam(3+3*i,pname,v_baselineshits[i].X(),0.02,-0.25,0.25);
mFit.SetParam(6+3*i,pname,v_baselineshits[i].X(),0.02,-0.25,0.25);
sprintf(pname,"BaselineShift_Y_%d",(int)(i+2));
mFit.SetParam(4+3*i,pname,v_baselineshits[i].Y(),0.02,-0.25,0.25);
mFit.SetParam(7+3*i,pname,v_baselineshits[i].Y(),0.02,-0.25,0.25);
sprintf(pname,"BaselineShift_Z_%d",(int)(i+2));
mFit.SetParam(5+3*i,pname,v_baselineshits[i].Z(),0.02,-0.25,0.25);
mFit.SetParam(8+3*i,pname,v_baselineshits[i].Z(),0.02,-0.25,0.25);
if (fgfixbaseline) {
cout << " ... fitting phases only, fixed baselines "<<endl;
mFit.SetFix(3+3*i); mFit.SetFix(4+3*i); mFit.SetFix(5+3*i);
mFit.SetFix(6+3*i); mFit.SetFix(7+3*i); mFit.SetFix(8+3*i);
}
}
cout << " Performing the fit (tot_npoints_fit= "<<tot_npoints_fit<<" ?= (npoints2="<<npoints2<<") ..."<< endl;
......@@ -894,14 +920,16 @@ int CxBaselineFitter::dofit(string outfilename, bool fgfixbaseline)
}
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_phi_0[i]=mFit.GetParm(2*i);
v_err_phi_0[i]=mFit.GetParmErr(2*i);
v_a_phi[i]=mFit.GetParm(2*i+1);
v_err_a_phi[i]=mFit.GetParmErr(2*i+1);
double xs=mFit.GetParm(i*3+6);
double exs=mFit.GetParmErr(i*3+6);
double ys=mFit.GetParm(i*3+7);
double eys=mFit.GetParmErr(i*3+7);
double zs=mFit.GetParm(i*3+8);
double ezs=mFit.GetParmErr(i*3+8);
v_baselineshits[i]=Vector3d(xs,ys,zs);
v_err_baselineshits[i]=Vector3d(exs,eys,ezs);
......@@ -920,7 +948,7 @@ 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);
My6CxMinZFunc mzfunc(v_acxd, v_trks, true);
MinZSimplex simplex(&mzfunc);
// Guess the center and step for constructing the initial simplex
size_t nparam = 4*(NB_ANTENNES-1);
......@@ -952,19 +980,19 @@ 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_phases[i]=oparm(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_baselineshits[i]=Vector3d(xs,ys,zs);
cout << " ANTENNE["<<i+2<<"] : Phase="<<v_phases[i]<<" BaseLineShift="<<v_baselineshits[i]<<endl;
cout << " ANTENNE["<<i+2<<"] : Phase="<<v_phi_0[i]<<" BaseLineShift="<<v_baselineshits[i]<<endl;
}
}
return 0;
}
int CxBaselineFitter::doCheck(bool fgminimize)
int CxBaselineFitter::doCheck()
{
size_t NB_ANTENNES=v_acxd[0].getNbAutoCor(); // nombre d'antennes
size_t NB_CXCORS=v_acxd[0].getNbCrossCor();
......@@ -974,7 +1002,7 @@ int CxBaselineFitter::doCheck(bool fgminimize)
if (NB_ANTENNES != 4)
throw PError("CxBaselineFitter::doCheck() NB_ANTENNES != 4 Current version works only for 4 antenna");
My6CxMinZFunc mzfunc(v_acxd, v_trks);
My6CxMinZFunc mzfunc(v_acxd, v_trks, true); // 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);
......
trkfit.h
View file @ 2511dec1
......@@ -42,6 +42,8 @@ public:
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
};
......@@ -80,7 +82,9 @@ public:
vector< vector<double> > v_min_cxdata;
vector< vector<double> > v_max_cxdata;