Commit 76f0a7b7 authored by Reza  ANSARI's avatar Reza ANSARI
Browse files

Fin debug programme de fit autocor et cross-cor des visibilites a partir de...

Fin debug programme de fit autocor et cross-cor des visibilites a partir de traces (satellites, ...) , Reza 22/12/2018
parent e78c0ea3
......@@ -24,7 +24,7 @@ public:
v_interp_elev(v_interp_elev_), v_interp_sazim(v_interp_sazim_),
nac_(nac), fggauss_beam_(fggauss_beam)
{
cout << " MyACSignal() Antenna=nac+1= "<<nac_+1<<" NPar="<< 3+2*v_time_data_.size() << endl;
//DBG 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++)
......@@ -33,9 +33,10 @@ public:
virtual int getDataSignal(size_t j, vector<double> & sig)
{
sig.resize(v_time_data[j].size());
//DEL 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];
//DEL for(size_t k=0; k<sig.size(); k++) sig[k]=vdata[k];
sig = vdata;
return (int)sig.size();
}
......@@ -43,8 +44,8 @@ public:
{
vector<double> sig;
getDataSignal(j, sig);
TVector<double> rvs(sig.size());
for(size_t k=0; k<sig.size(); k++) rvs(k)=sig[k];
TVector<double> rvs(sig);
//DEL for(size_t k=0; k<sig.size(); k++) rvs(k)=sig[k];
return rvs;
}
......
......@@ -21,14 +21,15 @@ public:
v_interp_elev(v_interp_elev_), v_interp_sazim(v_interp_sazim_),
ncx_(ncx), beam(beam_)
{
cout << " MyCxSignal() CxCorr=ncx "<<ncx_<< endl;
//DBG cout << " MyCxSignal() CxCorr=ncx "<<ncx_<< endl;
}
virtual int getDataSignal(size_t j, vector< complex<double> > & sig)
{
sig.resize(v_time_data[j].size());
//DEL sig.resize(v_time_data[j].size());
vector< complex<double> > & vcxdata = vv_cxdata[j][ncx_];
for(size_t k=0; k<sig.size(); k++) sig[k]=vcxdata[k];
sig = vcxdata;
//DEL for(size_t k=0; k<sig.size(); k++) sig[k]=vcxdata[k];
return (int)sig.size();
}
......@@ -36,8 +37,8 @@ public:
{
vector< complex<double> >sig;
getDataSignal(j, sig);
TVector< complex<double> > rvs(sig.size());
for(size_t k=0; k<sig.size(); k++) rvs(k)=sig[k];
TVector< complex<double> > rvs(sig);
//DEL for(size_t k=0; k<sig.size(); k++) rvs(k)=sig[k];
return rvs;
}
......@@ -57,7 +58,8 @@ public:
double phisrcdeg = 90.-sazim;
if (phisrcdeg < 0.) phisrcdeg+=360.;
complex<double> cxb=beam.Value(thetasrc, Angle(phisrcdeg,Angle::Degree).ToRadian());
sig[k] = complex<double>(A*cos(phase),sin(phase))*cxb;
// 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);
//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;
......@@ -96,7 +98,7 @@ public:
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];
xx=(sig[k].imag()-vcxdata[k].imag())/verr[k];
xi2 += (xx*xx);
}
}
......
......@@ -2,7 +2,8 @@
PAON4 analysis :
Determining antenna pointing and dish size fitting parameters ...
(source track corresponds to a satellite track (or planet or sun ...)
Uses gacfit.h and acfit.h
Does also fit of phases for the 6 cross-correlations
Uses gacfit.h and gcxfit.h acbeam.h cxbeam.h
R. Ansari , December 2018
Example command to run :
......@@ -56,11 +57,17 @@ static bool fggaussbeam=true;
// --- if true, perform Beam normalisation
// static bool fgdobeamnormalise=false;
static string outfilename; // output filename for fit parameters
static string outcheckfilename; // output filename for fit parameters
static string outfilename; // output filename for auto-correlation fit parameters
static string outfilenamecx; // output filename for cross-correlation fit parameters
static string outcheckfilename; // output filename for PAON4 and computed auto-cor using fitted parameters
static string outcheckfilenamecx; // output filename for PAON4 and computed cross-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=1; // print level
static int prtlevel=0; // print level
//---- input data definition
static size_t NB_ANTENNES = 4; // Number of antenna ( 4 for PAON4 )
......@@ -102,7 +109,8 @@ int do_p4_cxtrkfit(vector< vector<double> > & v_time_data, vector< vector< vecto
vector< vector<double> > & v_min_cxdata, vector< vector<double> > & v_max_cxdata,
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);
vector<double> & v_phase, vector< vector<double> > & v_Acx,
vector< vector<double> > & v_A, bool useAac=true);
//--------------------------------------------------------------
//--------- Main program -------------
......@@ -146,36 +154,10 @@ int main (int narg, char* arg[])
int rcf = do_p4_trkfit(v_time_data, vv_data, vv_err, v_freqs, v_min_data, v_max_data,
v_interp_elev, v_interp_sazim, v_Ddish, v_thetaant, v_phiant, v_A, v_B);
if (outcheckfilename.length()>0) {
cout << " 5.a Check : computing expected signal for fitted params , will be saved to file "<<outcheckfilename<<endl;
POutPersist pos(outcheckfilename);
for(size_t ii=0; ii<4; ii++) {
cout << " 5.b Computing DataSignal & Expected Signal for fitted params and dish "<<ii+1<<endl;
MyACSignal macs(v_time_data, vv_data, vv_err, v_freqs, v_interp_elev, v_interp_sazim, ii, fggaussbeam);
double Ddishfit=v_Ddish[ii];
double thetafit=v_thetaant[ii];
double phifit=v_phiant[ii];
char oname[32];
for(size_t j=0; j<NTRK; j++) {
double A = v_A[ii][j];
double B = v_B[ii][j];
Vector signal = macs.getDataSignal(j);
sprintf(oname,"data_%d_%d",(int)ii+1,(int)j+1);
pos << PPFNameTag(oname)<<signal;
Vector expsignal = macs.getExpectedSignal(j, Ddishfit, thetafit, phifit, A, B);
sprintf(oname,"simu_%d_%d",(int)ii+1,(int)j+1);
pos << PPFNameTag(oname)<<expsignal;
}
}
} //End-of-check , saving data & expected signal to outcheckfilename
vector<double> v_phase(NB_CXCORS);
vector< vector<double> > v_Acx(NB_CXCORS);
int rcf2 = do_p4_cxtrkfit(v_time_data, vv_cxdata, vv_cxerr, v_freqs, v_min_cxdata, v_max_cxdata,
v_interp_elev, v_interp_sazim, v_Ddish, v_thetaant, v_phiant, v_phase, v_Acx);
v_interp_elev, v_interp_sazim, v_Ddish, v_thetaant, v_phiant, v_phase, v_Acx, v_A, fguseAac4Cx);
} // End of try bloc
catch (PThrowable & exc) {
......@@ -265,7 +247,7 @@ int read_p4_autocxcors(vector< vector<double> > & v_time_data, vector< vector< v
vector< complex<double> > & vcx = v_vcx[ii];
v_cxdata[ii].push_back(vcx[k]);
double acx=std::abs(vcx[k]);
v_cxerr[ii].push_back(0.1*acx);
v_cxerr[ii].push_back(0.1*sqrt(acx));
if (acx<v_min_cxdata[j][ii]) v_min_cxdata[j][ii]=acx;
if (acx>v_max_cxdata[j][ii]) v_max_cxdata[j][ii]=acx;
}
......@@ -346,8 +328,7 @@ 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)
{
cout << "--------------------------------------------------------------------------------------"<<endl;
cout << "--------------------------------------------------------------------------------------"<<endl;
cout << "======================================================================================"<<endl;
cout << "---- trkacfit/do_p4_trkfit() ; Performing antenna pointing fit ..."<<endl;
ofstream ofr(outfilename.c_str());
ofr << "#### Pointing/dish diameter fit on autocorrelation (programme trkacfit.cc ) "<<endl
......@@ -358,7 +339,7 @@ int do_p4_trkfit(vector< vector<double> > & v_time_data, vector< vector< vector<
int tot_npoints = 0;
for(size_t j=0; j<NTRK; j++) tot_npoints += v_time_data[j].size();
for(size_t ii=0; ii<NB_ANTENNES; ii++) {
cout << "do_p4_trkfit 1."<<ii+1<<" Creating General Fit for AutoCor Antenna= " << ii+1 << endl;
cout << "-------- do_p4_trkfit 1."<<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];
......@@ -403,12 +384,12 @@ int do_p4_trkfit(vector< vector<double> > & v_time_data, vector< vector< vector<
}
//DBG mFit.PrintFit();
cout << "do_p4_trkfit 2."<<ii+1<<" Performing the fit for AutoCor Antenna= " << ii+1 << endl;
// cout << "do_p4_trkfit 2."<<ii+1<<" Performing the fit for AutoCor Antenna= " << ii+1 << endl;
int rcfit = mFit.Fit();
if (prtlevel>1) mFit.PrintFit();
if(rcfit>0) {
cout<< "-------------------------- Result for Antenna No " << ii+1 << endl;
cout<< "---Fit resultReduce_Chisquare = " << mFit.GetChi2Red()<< " nstep="<<mFit.GetNStep() << " rc="<<rcfit<<endl;
cout<< "------- Fit result for Antenna No="<<ii+1<<" Reduce_Chisquare = " << mFit.GetChi2Red()
<< " nstep="<<mFit.GetNStep() << " rc="<<rcfit<<endl;
ofr <<setw(4)<<ii+1<<" "<<setw(8)<<mFit.GetChi2Red()<<" ";
double Dfit=mFit.GetParm(0); double err_Dfit=mFit.GetParmErr(0);
cout <<setw(16)<<"DishDiameter= "<<setw(10)<<Dfit<<" +/- "<<setw(10)<<err_Dfit<<" m."<<endl;
......@@ -441,12 +422,38 @@ int do_p4_trkfit(vector< vector<double> > & v_time_data, vector< vector< vector<
ofr << endl;
}
else {
cout << "--- Fit_Error, rc = " << rcfit << " nstep="<<mFit.GetNStep()<<endl;
cout << "---Fit failed for "<<ii+1<<"--- 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);
}
}
if (outcheckfilename.length()>0) {
cout << "-----do_p4_trkfit 3/ Check : computing expected signal for fitted params , will be saved to file "<<outcheckfilename<<endl;
POutPersist pos(outcheckfilename);
for(size_t ii=0; ii<4; ii++) {
if (prtlevel>1) cout << "do_p4_trkfit 3.b/ Computing DataSignal & Expected Signal for fitted params and dish "<<ii+1<<endl;
MyACSignal macs(v_time_data, vv_data, vv_err, v_freqs, v_interp_elev, v_interp_sazim, ii, fggaussbeam);
double Ddishfit=v_Ddish[ii];
double thetafit=v_thetaant[ii];
double phifit=v_phiant[ii];
char oname[32];
for(size_t j=0; j<NTRK; j++) {
double A = v_A[ii][j];
double B = v_B[ii][j];
Vector signal = macs.getDataSignal(j);
sprintf(oname,"ac_%d_%d",(int)ii+1,(int)j+1);
pos << PPFNameTag(oname)<<signal;
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;
}
}
} //End-of-check , saving data & expected signal to outcheckfilename
return 0;
}
......@@ -456,25 +463,30 @@ int do_p4_cxtrkfit(vector< vector<double> > & v_time_data, vector< vector< vecto
vector< vector<double> > & v_min_cxdata, vector< vector<double> > & v_max_cxdata,
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 << "--------------------------------------------------------------------------------------"<<endl;
cout << "--------------------------------------------------------------------------------------"<<endl;
vector<double> & v_phase, vector< vector<double> > & v_Acx,
vector< vector<double> > & v_A, bool useAac)
cout << "---- trkacfit/do_p4_cxtrkfit() ; Performing cross-cor phase fit ..."<<endl;
POutPersist pox("cxdbg.ppf");
ofstream ofr("cxcorfit.txt");
{
cout << "======================================================================================"<<endl;
cout << "---------- trkacfit/do_p4_cxtrkfit() ; Performing cross-cor phase fit ..."<<endl;
if (useAac) cout << " ... Using Amplitude from auto-correlations fit for initial fit parameter value..."<<endl;
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);
ofr << "#### cross-cor phase fit (programme trkacfit.cc ) "<<endl
<< "## NumCxCor Xi2red Phase err_Phase A0 err_A0 A1 err_A1 ..."<<endl;
<< "## NumCxCor Xi2red Phase err_Phase (deg) 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++) {
v_Acx[ii].resize(NTRK);
Vector3d baseline=P4Coords::getBaseline(Anum1[ii]+1,Anum2[ii]+1);
cout << "do_p4_cxtrkfit() 1."<<ii+1<<" Creating General Fit for CrossCor= " << ii << " FxF= " << Anum1[ii]+1<<"x"<<Anum2[ii]+1
<<" Baseline="<<baseline<<endl;
cout << "--------- 1."<<ii+1<<" do_p4_cxtrkfit() Doing fit for CrossCor= " << ii << " FxF= "
<< Anum1[ii]+1<<"x"<<Anum2[ii]+1<<" Baseline="<<baseline<<endl;
GeneralFitData gdata(1, tot_npoints);
for(size_t j=0; j<NTRK; j++) {
vector< vector< complex<double> > > & v_cxdata = vv_cxdata[j];
......@@ -492,11 +504,6 @@ int do_p4_cxtrkfit(vector< vector<double> > & v_time_data, vector< vector< vecto
acb2.setGaussianLobe(fggaussbeam);
CxBeam cxbeam(acb1, acb2, baseline);
/*
cout << "do_p4_cxtrkfit()*DBG*A0* Theta_0= "<<theta_0<<" Phi_0= "<<phi_0<<endl;
cout << "do_p4_cxtrkfit()*DBG*A1* Theta1="<<v_thetaant[Anum1[ii]]<<" Phi1="<<v_phiant[Anum1[ii]]<< " Beam1="<<acb1.Value(theta_0, phi_0)<<endl;
cout << "do_p4_cxtrkfit()*DBG*A2* Theta2="<<v_thetaant[Anum2[ii]]<<" Phi2="<<v_phiant[Anum2[ii]]<< " Beam2="<<acb2.Value(theta_0, phi_0)<<endl;
*/
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
......@@ -512,18 +519,21 @@ int do_p4_cxtrkfit(vector< vector<double> > & v_time_data, vector< vector< vecto
Vector asig = SOPHYA::abs(signal);
double mins, maxs;
asig.MinMax(mins, maxs);
sprintf(oname,"cx_%d_%d",(int)ii+1,(int)j+1);
pox << PPFNameTag(oname)<<signal;
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;
aexpsig.MinMax(mine, maxe);
A=maxs/maxe;
v_amp[j]=A;
expsignal *= complex<double>(A, 0.);
//DBG cout << "do_p4_cxtrkfit()*DBG*B* j="<<j<<" A="<<A<<" max-sig= "<<maxs<<" max-exp="<<maxe<<endl;
sprintf(oname,"simcx_%d_%d",(int)ii+1,(int)j+1);
pox << PPFNameTag(oname)<<expsignal;
/*DBG if (pox) {
expsignal *= complex<double>(A,0.);
sprintf(oname,"ecx_%d_%d",(int)ii+1,(int)j+1);
(*pox) << PPFNameTag(oname)<<expsignal;
} */
}
double fparm[50]; fparm[0]=0.;
......@@ -532,10 +542,12 @@ int do_p4_cxtrkfit(vector< vector<double> > & v_time_data, vector< vector< vecto
double bestphase=0.;
int bestnpts,npts;
int bestafact;
double afact[12]={0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.35,1.5,1.75,2.0};
double afact[10]={0.3,0.5,0.75,1.0,1.25,1.5,1.75,2.0,2.5,3.};
for(int ia=0; ia<12; ia++) {
for(size_t j=0; j<NTRK; j++) {
fparm[1+3*j]=afact[ia]*v_amp[j]; fparm[2+3*j]=fparm[3+3*j]=0.;
double Aac=sqrt(v_A[Anum1[ii]][j] * v_A[Anum2[ii]][j]);
fparm[1+3*j]=(useAac?Aac:v_amp[j]);
fparm[1+3*j]*=afact[ia]; fparm[2+3*j]=fparm[3+3*j]=0.;
}
for(double ph=0.; ph<360.; ph += 1) {
fparm[0]=Angle(ph, Angle::Degree).ToRadian();
......@@ -545,16 +557,23 @@ int do_p4_cxtrkfit(vector< vector<double> > & v_time_data, vector< vector< vecto
}
}
}
mFit.SetParam(0,"Phase",bestphase,M_PI/360.,0.,2.2*M_PI);
cout << "do_p4_cxtrkfit() 2."<<ii+1<<" bestxi2="<<bestxi2<<" bestphase="<<Angle(bestphase).ToDegree()<<" bestnpts="<<bestnpts
<<" bestafact="<<bestafact<< " A= "; for(size_t j=0; j<NTRK; j++) cout << v_amp[j] << " , "; cout << endl;
mFit.SetParam(0,"Phase",bestphase,M_PI/720.,0.,2.2*M_PI);
cout << "2."<<ii+1<<" bestxi2="<<bestxi2<<" bestphase="<<Angle(bestphase).ToDegree()
<<" bestnpts="<<bestnpts<<" bestafact="<<bestafact<< " A= ";
v_phase[ii]=bestphase;
for(size_t j=0; j<NTRK; j++) {
cout << v_amp[j] << " , "; cout << endl;
double Aac=sqrt(v_A[j][Anum1[ii]] * v_A[j][Anum2[ii]]);
v_Acx[ii][j]=v_amp[j];
}
for(size_t j=0; j<NTRK; j++) {
char pname[32];
sprintf(pname,"A%d",(int)(j+1));
double A=v_amp[j]; // v_max_cxdata[j][ii];
// double A = sqrt(v_A[j][Anum1[ii]] * v_A[j][Anum2[ii]]);
mFit.SetParam(1+3*j,pname,A,A/10.,A/4,A*3);
sprintf(pname,"A%d",(int)(j+1));
double Aac=sqrt(v_A[Anum1[ii]][j] * v_A[Anum2[ii]][j]);
double A=(useAac?Aac:v_amp[j]);
//DBG cout << "*DBG* j="<<j<<" Aac= "<<Aac<<" v_amp="<<v_amp[j]<<" A= "<<A<<" A1="<<v_A[Anum1[ii]][j]<<" A2="<<v_A[Anum2[ii]][j]<<endl;
mFit.SetParam(1+3*j,pname,A,A/10.,A/4,A*4);
sprintf(pname,"Bre%d",(int)(j+1));
mFit.SetParam(2+3*j,pname,0.,A/25.,-A/5,A/5.);
sprintf(pname,"Bim%d",(int)(j+1));
......@@ -563,32 +582,43 @@ int do_p4_cxtrkfit(vector< vector<double> > & v_time_data, vector< vector< vecto
mFit.SetFix(3+3*j,0.);
}
//DBG mFit.PrintFit();
cout << "do_p4_cxtrkfit() 3."<<ii+1<<" Performing the fit for CrossCor " << ii << " FxF= " << Anum1[ii]+1<<"x"<<Anum2[ii]+1<<endl;
if (prtlevel>1)
cout << " 3."<<ii+1<<" Performing the fit for CrossCor " << ii << " FxF= " << Anum1[ii]+1<<"x"<<Anum2[ii]+1<<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;
// cout<< "-------------------------- Result for Cross No " << ii << endl;
cout<< "------ Fit result for Cross No "<<ii<<" Reduce_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)<<Angle(phase).ToDegree()<<" +/- "<<setw(10)<<Angle(err_phase).ToDegree()<<" deg."<<endl;
ofr <<setw(8)<<phase<<" "<<setw(8)<<err_phase<<" ";
// v_Ddish[ii]=Dfit;
cout <<"Phase= "<<setw(10)<<Angle(phase).ToDegree()<<" +/- "<<setw(10)<<Angle(err_phase).ToDegree()<<" deg."<<endl;
ofr <<setw(8)<<Angle(phase).ToDegree()<<" "<<setw(8)<<Angle(err_phase).ToDegree()<<" ";
v_phase[ii]=phase;
for(size_t j=0; j<NTRK; j++) {
double Aac=sqrt(v_A[Anum1[ii]][j] * v_A[Anum2[ii]][j]);
double Ai=(useAac?Aac:v_amp[j]);
double A=mFit.GetParm(1+3*j); double err_A=mFit.GetParmErr(1+3*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<<" ";
cout << " Trk["<<j<<"] A= "<<A<<" +/- "<<err_A<<" (A/Ai="<<A/Ai<<")"<<endl;
v_Acx[ii][j]=A; // v_B[ii][j]=B;
ofr <<setw(8)<<A<<" "<<setw(8)<<err_A<<" ";
}
ofr << endl;
}
else {
cout << "--- Fit_Error, rc = " << rcfit << " nstep="<<mFit.GetNStep()<<endl;
cout << "---Fit failed for "<<ii<<" 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);
}
}
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;
return 0;
}
......@@ -628,7 +658,9 @@ int decode_args(int narg, char* arg[])
}
outfilename = "trkacfit.txt";
outfilenamecx = "trkacfit_cx.txt";
outcheckfilename = "";
outcheckfilenamecx = "";
input_base_path = "./";
prtlevel=0;
D_dish = 4.5;
......@@ -673,18 +705,28 @@ int decode_args(int narg, char* arg[])
return 9;
}
size_t pp=outfilename.rfind('.');
if (pp!=string::npos) outfilenamecx=outfilename.substr(0,pp)+"_cx"+outfilename.substr(pp);
else outfilenamecx=outfilename+"_cx";
NTRK=lastargs.size();
cout << " ---------- trkacfit/run parameters:"<<endl;
cout << " Beam: D_dish (initial guess)="<<D_dish<<" GaussianBeam ? "<<(fggaussbeam?"Yes":"No")
<< " dec-shift(zenithAngle)= "<<zenang<<endl;
cout << " Input Base Path= "<< input_base_path << " NTrack="<<NTRK<<" OutFileName= "<<outfilename<<endl;
cout << " Input Base Path= "<< input_base_path << " NTrack="<<NTRK<<" OutFileName= "<<outfilename
<< " OutFileNameCx= "<<outfilenamecx<<endl;
if (outcheckfilename.length()>0) {
pp=outcheckfilename.rfind('.');
if (pp!=string::npos) outcheckfilenamecx=outcheckfilename.substr(0,pp)+"_cx"+outcheckfilename.substr(pp);
else outcheckfilenamecx=outcheckfilename+"_cx";
cout << " data and expected signals with be saved to files "<<outcheckfilename<<" , Cx:"<<outcheckfilenamecx<<endl;
}
if (zenang<0.) {
theta_0 = Angle(-zenang, Angle::Degree).ToRadian(); phi_0 = Angle::PioTwoCst()+Angle::OnePiCst();
}
else {
theta_0 = Angle(+zenang, Angle::Degree).ToRadian(); phi_0 = Angle::PioTwoCst();
}
NTRK=lastargs.size();
dataflnm.resize(NTRK);
tstart.resize(NTRK); tend.resize(NTRK);
v_freqs.resize(NTRK);
......@@ -710,7 +752,8 @@ int decode_args(int narg, char* arg[])
v_freqs[i]=freq;
cout << "--DataSet["<<i+1<<"] Data="<<dataflnm[i]<<" Trk= "<<trkflnm[i]<<" "<<tstart[i]<<" <=t<= "<<tend[i]
<< "("<<t_s<<" <=t_min<= "<<t_e<<" @ Freq= "<<v_freqs[i]<<" MHz"<<endl;
}
}
cout<<"---------------------------------------------------------------------------------"<<endl;
return 0;
}
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