Commit 0b361dd2 authored by Reza  ANSARI's avatar Reza ANSARI
Browse files

Suite implementation fit de phase sur les cross-correlations ds trkacfit.cc , Reza 21/12/2018

parent 824c1e6f
......@@ -121,13 +121,15 @@ public:
// circular beam response
double alp=acos(bdir_.Psc(uvo))*DoL_;
if (fggauss_) {
return ( exp(-2.*alp*alp)*NormFac_ ); // Jiao : Check the factor 2 in the exp( )
// JEC , wiki BAORadio 2.1767 = 0.22054*Pi^2
return ( exp(-2.1767*alp*alp)*NormFac_ ); // Jiao : Check the factor 2 in the exp( ) -> exp(-2.*alp*alp)
}
else {
if (alp<1.e-19) return NormFac_;
alp*=M_PI; double xa=2.*j1(alp)/alp;
return (NormFac_*xa*xa);
}
}
// Return the beam in the form of HEALPix spherical map
......
......@@ -26,7 +26,7 @@ public:
{
}
// copy operator
CxBeam& operator = (CxBeam const& a)
CxBeam& operator = (CxBeam const& c)
{
a1_=c.a1_; a2_=c.a2_; baseline_=c.baseline_; facphase_=c.facphase_;
return (*this);
......@@ -53,3 +53,5 @@ public:
Vector3d baseline_;
double facphase_;
};
#endif
......@@ -14,45 +14,41 @@
class MyCxSignal {
public:
MyCxSignal(vector< vector<double> > & v_time_data_, vector< vector< vector< complex<double> > > > & vv_cxdata_,
vector< vector< vector<double> > > & vv_cxerr_, vector<double> & v_freqs_,
vector< vector< vector<double> > > & vv_cxerr_, // vector<double> & v_freqs_,
vector< SLinInterp1D > & v_interp_elev_, vector< SLinInterp1D > & v_interp_sazim_,
size_t ncx=0, bool fggauss_beam=true)
: v_time_data(v_time_data_), vv_data(vv_data_), vv_err(vv_err_), v_freqs(v_freqs_),
CxBeam & beam_, size_t ncx=0)
: 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_),
nac_(nac),
ncx_(ncx), beam(beam_)
{
cout << " MyACSignal() Antenna=nac+1= "<<nac_+1<<" NPar="<< 3+2*v_time_data_.size() << endl;
double clight = PhysQty::c().SIValue();
v_lambdas.resize(v_freqs.size());
for(size_t k=0; k<v_freqs.size(); k++)
v_lambdas[k] = clight/(v_freqs[k]*1.e6);
cout << " MyCxSignal() CxCorr=ncx "<<ncx_<< endl;
}
virtual int getDataSignal(size_t j, vector<double> & sig)
virtual int getDataSignal(size_t j, vector< complex<double> > & sig)
{
sig.resize(v_time_data[j].size());
vector<double> & vdata = vv_data[j][nac_];
for(size_t k=0; k<sig.size(); k++) sig[k]=vdata[k];
vector< complex<double> > & vcxdata = vv_cxdata[j][ncx_];
for(size_t k=0; k<sig.size(); k++) sig[k]=vcxdata[k];
return (int)sig.size();
}
virtual TVector<double> getDataSignal(size_t j)
{
vector<double> sig;
vector< complex<double> >sig;
getDataSignal(j, sig);
TVector<double> rvs(sig.size());
TVector< complex<double> > rvs(sig.size());
for(size_t k=0; k<sig.size(); k++) rvs(k)=sig[k];
return rvs;
}
virtual int getExpectedSignal(size_t j, vector<double> & sig, double D_dish, double theta_b, double phi_b, double A, double B)
virtual int getExpectedSignal(size_t j, vector< complex<double> > & sig, double phase, double A)
{
//DBG if (debug_gacfit > 0) cout<<"*DBG*A*getExpectedSignal() j="<<j<<" size="<<v_time_data[j].size()<<endl;
sig.resize(v_time_data[j].size());
vector<double> & vtime = v_time_data[j];
vector<double> & vdata = vv_data[j][nac_];
vector<double> & verr = vv_err[j][nac_];
// vector< complex<double> > & vcxdata = vv_cxdata[j][nac_];
// vector<double> & vcxerr = vv_cxerr[j][nac_];
for(size_t k=0; k<vtime.size(); k++) {
double elev=v_interp_elev[j].YInterp(vtime[k]);
double thetasrc=Angle(90.-elev,Angle::Degree).ToRadian();
......@@ -60,12 +56,8 @@ public:
while (sazim > 360.) sazim -= 360.;
double phisrcdeg = 90.-sazim;
if (phisrcdeg < 0.) phisrcdeg+=360.;
// double phisrc=acos(v_interp_cphi[j].YInterp(vtime[k]));
// if (phisrc<0.) phisrc+=Angle::TwoPiCst();
ACBeam beam(D_dish, theta_b, phi_b, v_lambdas[j]);
beam.setGaussianLobe(fggauss_beam_);
double acb=beam.Value(thetasrc, Angle(phisrcdeg,Angle::Degree).ToRadian());
sig[k] = A*acb+B;
complex<double> cxb=beam.Value(thetasrc, Angle(phisrcdeg,Angle::Degree).ToRadian());
sig[k] = complex<double>(A*cos(phase),sin(phase))*cxb;
//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;
......@@ -73,30 +65,20 @@ public:
return (int)sig.size();
}
virtual int getExpectedSignal(size_t j, vector<double> & sig, const double* parm)
virtual int getExpectedSignal(size_t j, vector< complex<double> > & sig, const double* parm)
{
/*
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
*/
double D_dish = parm[0];
double theta_b = parm[1];
double phi_b = parm[2];
double A = parm[3+2*j];
double B = parm[4+2*j];
return getExpectedSignal(j, sig, D_dish, theta_b, phi_b, A, B);
double phase = parm[0];
double A = parm[1+j];
return getExpectedSignal(j, sig, phase, A);
}
virtual TVector<double> getExpectedSignal(size_t j, double D_dish, double theta_b, double phi_b, double A, double B)
virtual TVector< complex<double> > getExpectedSignal(size_t j, double phase, double A)
{
vector<double> sig;
vector< complex<double> > sig;
//DBG cout << "*DBG*getExpectedSignal() : calling getExpectedSignal(j="<<j<<" ,sig,...)"<<endl;
//DBG debug_gacfit=1;
getExpectedSignal(j, sig, D_dish, theta_b, phi_b, A, B);
TVector<double> rvs(sig.size());
getExpectedSignal(j, sig, phase, A);
TVector< complex<double> > rvs(sig.size());
for(size_t k=0; k<sig.size(); k++) rvs(k)=sig[k];
return rvs;
}
......@@ -106,25 +88,51 @@ public:
double xi2=0.;
npts=0;
for(size_t j=0; j<v_time_data.size(); j++) {
vector<double> sig;
vector< complex<double> > sig;
npts += getExpectedSignal(j, sig, parm);
vector<double> & vdata = vv_data[j][nac_];
vector<double> & verr = vv_err[j][nac_];
for(size_t k=0; k<vdata.size(); k++) {
double xx=(sig[k]-vdata[k])/verr[k];
vector< complex<double> > & vcxdata = vv_cxdata[j][ncx_];
vector<double> & verr = vv_cxerr[j][ncx_];
for(size_t k=0; k<vcxdata.size(); k++) {
double xx=(sig[k].real()-vcxdata[k].real())/verr[k];
xi2 += (xx*xx);
xx=(-sig[k].imag()-vcxdata[k].imag())/verr[k];
xi2 += (xx*xx);
}
}
npts *= 2;
return xi2;
}
vector< vector<double> > & v_time_data;
vector< vector< vector<double> > > & vv_data;
vector< vector< vector<double> > > & vv_err;
vector<double> & v_freqs;
vector<double> v_lambdas;
vector< vector< vector< complex<double> > > > & vv_cxdata;
vector< vector< vector<double> > > & vv_cxerr;
vector< SLinInterp1D > & v_interp_elev;
vector< SLinInterp1D > & v_interp_sazim;
size_t nac_;
bool fggauss_beam_;
CxBeam beam;
size_t ncx_;
};
//-----------------------------------------------------------------
class MyCxGenXi2 : public GeneralXi2, MyCxSignal {
public:
MyCxGenXi2(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_,
CxBeam & beam, size_t ncx=0)
: GeneralXi2(1+v_time_data_.size()) ,
MyCxSignal(v_time_data_, vv_cxdata_, vv_cxerr_, v_interp_elev_, v_interp_sazim_, beam, ncx)
{
}
virtual double Value(GeneralFitData& data, double* parm, int& ndataused)
/*
param[0] = phase
parm[1+j] = Aj (Amplitude) for track j
*/
{
return getXi2(parm, ndataused);
}
};
#endif
......@@ -90,7 +90,7 @@ trkacfit : $(EXE)/trkacfit
$(EXE)/trkacfit : $(OBJ)/trkacfit.o $(MYOLISTHERE)
$(CXXLINK) -o $(EXE)/trkacfit $(OBJ)/trkacfit.o $(MYOLISTHERE) $(SOPHYAEXTSLBLIST)
$(OBJ)/trkacfit.o : trkacfit.cc $(MYINCLISTHERE) acbeam.h gacfit.h
$(OBJ)/trkacfit.o : trkacfit.cc $(MYINCLISTHERE) acbeam.h gacfit.h cxbeam.h gcxfit.h
$(CXXCOMPILE) -o $(OBJ)/trkacfit.o trkacfit.cc
## programme pour extraire un DataTable a partir des TFM (fichier PPF), input programme trkacfit
......
......@@ -32,8 +32,11 @@
#include "fitsioserver.h"
#include "slininterp.h"
#include "p4autils.h"
#include "acbeam.h"
#include "gacfit.h"
#include "gcxfit.h"
using namespace std;
......@@ -92,6 +95,12 @@ int do_p4_trkfit(vector< vector<double> > & v_time_data, vector< vector< vector<
vector<double> & v_Ddish, vector<double> & v_thetaant, vector<double> & v_phiant,
vector< vector<double> > & v_A, vector< vector<double> > & v_B);
int do_p4_cxtrkfit(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<double> & v_Ddish, vector<double> & v_thetaant, vector<double> & v_phiant,
vector<double> & v_phase, vector< vector<double> > & v_A);
//--------------------------------------------------------------
//--------- Main program -------------
//--------------------------------------------------------------
......@@ -158,6 +167,10 @@ int main (int narg, char* arg[])
}
} //End-of-check , saving data & expected signal to outcheckfilename
vector<double> v_phase;
int rcf2 = do_p4_cxtrkfit(v_time_data, vv_cxdata, vv_cxerr, v_freqs,
v_interp_elev, v_interp_sazim, v_Ddish, v_thetaant, v_phiant, v_phase, v_A);
} // End of try bloc
catch (PThrowable & exc) {
cerr << " trkacfit.cc: Catched Exception (PThrowable)" << (string)typeid(exc).name()
......@@ -419,6 +432,82 @@ int do_p4_trkfit(vector< vector<double> > & v_time_data, vector< vector< vector<
return 0;
}
/* --Fonction-- */
int do_p4_cxtrkfit(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<double> & v_Ddish, vector<double> & v_thetaant, vector<double> & v_phiant,
vector<double> & v_phase, vector< vector<double> > & v_A)
{
cout << "---- trkacfit/do_p4_cxtrkfit() ; Performing cross-cor phase fit ..."<<endl;
ofstream ofr("cxcorfit.txt");
ofr << "#### cross-cor phase fit (programme trkacfit.cc ) "<<endl
<< "## NumCxCor Xi2red Phase err_Phase A0 err_A0 A1 err_A1 ..."<<endl;
int tot_npoints = 0;
for(size_t j=0; j<NTRK; j++) tot_npoints += 2*(v_time_data[j].size());
size_t Anum1[6]={0,0,0,1,1,2};
size_t Anum2[6]={1,2,3,2,3,3};
for(size_t ii=0; ii<NB_CXCORS; ii++) {
cout << "3."<<ii+1<<" Creating General Fit for AutoCor Antenna= " << ii+1 << endl;
GeneralFitData gdata(1, tot_npoints);
/* for(size_t j=0; j<NTRK; j++) {
vector< vector<double> > & v_data = vv_data[j];
vector< vector<double> > & v_err = vv_err[j];
for(size_t k=0; k<v_time_data[j].size(); k++) {
gdata.AddData1(v_time_data[j][k],v_data[ii][k],v_err[ii][k]); // Fill x, y and error on y
}
} */
double clight = PhysQty::c().SIValue();
double lambda = clight/(v_freqs[0]*1.e6);
ACBeam acb1(v_Ddish[Anum1[ii]], v_thetaant[Anum1[ii]], v_phiant[Anum1[ii]], lambda);
ACBeam acb2(v_Ddish[Anum2[ii]], v_thetaant[Anum2[ii]], v_phiant[Anum2[ii]], lambda);
Vector3d baseline=P4Coords::getBaseline(Anum1[ii]+1,Anum2[ii]+1);
CxBeam cxbeam(acb1, acb2, baseline);
MyCxGenXi2 gxi2( v_time_data, vv_cxdata, vv_cxerr, v_interp_elev, v_interp_sazim, cxbeam, ii);
GeneralFit mFit(&gxi2);
mFit.SetData(&gdata); // connect data to the fitter , here the data is unused - gxi2 includes its data
mFit.SetMaxStep(1000);
// SetParam(int n,double value, double step,double min=1., double max=-1.);
mFit.SetParam(0,"Phase",M_PI,M_PI/360.,0.,2.2*M_PI);
for(size_t j=0; j<NTRK; j++) {
char pname[32];
sprintf(pname,"A%d",(int)(j+1));
double A = sqrt(v_A[j][Anum1[ii]] * v_A[j][Anum2[ii]]);;
mFit.SetParam(1+j,pname,A,A/10.,A/2,A*2);
}
//DBG mFit.PrintFit();
cout << "4."<<ii+1<<" Performing the fit for CrossCor " << ii << endl;
int rcfit = mFit.Fit();
if (prtlevel>1) mFit.PrintFit();
if(rcfit>0) {
cout<< "-------------------------- Result for Cross No " << ii << endl;
cout<< "---Fit resultReduce_Chisquare = " << mFit.GetChi2Red()<< " nstep="<<mFit.GetNStep() << " rc="<<rcfit<<endl;
ofr <<setw(4)<<ii+1<<" "<<setw(8)<<mFit.GetChi2Red()<<" ";
double phase=mFit.GetParm(0); double err_phase=mFit.GetParmErr(0);
cout <<setw(16)<<"Phase= "<<setw(10)<<phase<<" +/- "<<setw(10)<<err_phase<<" m."<<endl;
ofr <<setw(8)<<phase<<" "<<setw(8)<<err_phase<<" ";
// v_Ddish[ii]=Dfit;
for(size_t j=0; j<NTRK; j++) {
double A=mFit.GetParm(1+j); double err_A=mFit.GetParmErr(1+j);
cout << " Trk/Sat["<<j<<"] -> A= "<<A<<" +/- "<<err_A<<endl;
// v_A[ii][j]=A; v_B[ii][j]=B;
// ofr <<setw(8)<<A<<" "<<setw(8)<<err_A<<" "<<setw(8)<<B<<" "<<setw(8)<<err_B<<" ";
}
ofr << endl;
}
else {
cout << "--- Fit_Error, rc = " << rcfit << " nstep="<<mFit.GetNStep()<<endl;
ofr <<setw(4)<<ii+1<<" ERROR FIT RC="<<rcfit<<" nstep="<<mFit.GetNStep()<<endl;
if (prtlevel>0) mFit.PrintFitErr(rcfit);
}
}
return 0;
}
/* --Fonction-- */
int decode_args(int narg, char* arg[])
//! Decode program arguments
......
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