Commit 63b47b0a authored by Reza  ANSARI's avatar Reza ANSARI
Browse files

Implementation de toutes de les fonctions de lecture des donnees et de fit de...

Implementation de toutes de les fonctions de lecture des donnees et de fit de trkacfit.cc ds trkfit.h .cc , Reza 16/02/2019
parent ae86d4b8
......@@ -4,10 +4,28 @@
R. Ansari, Fevrier 2019 */
#include <iomanip>
#include "pexceptions.h"
#include "trkfit.h"
#include "datacards.h"
#include "skymap.h"
#include "array.h"
#include "acbeam.h"
#include "gacfit.h"
#include "gcxfit.h"
#include "gcxfitbaseline.h"
#include "p4autils.h"
//------------------- Print Level for this file --------------------------
static int _prtlevel_ =0;
void TrkFit_SetPrintLevel(int lev)
{
_prtlevel_=lev;
return;
}
//------------------- TrkInputDataSet -------------------------------------
......@@ -110,3 +128,536 @@ ostream & TrkInputDataSet::Print(ostream & os) const
}
//------------------------ ACxDataSet -------------------------------------
AcxDataSet::AcxDataSet(TrkInputDataSet & tkds)
: tot_npoints(0),zenang(0.),theta_0(0.),phi_0(0.)
{
ReadData(tkds);
}
AcxDataSet::AcxDataSet(AcxDataSet const & a)
: v_time_data(a.v_time_data), vv_data(a.vv_data), vv_err(a.vv_err),
v_min_data(a.v_min_data), v_max_data(a.v_max_data),
vv_cxdata(a.vv_cxdata), v_min_cxdata(a.v_min_cxdata), v_max_cxdata(a.v_max_cxdata),
tot_npoints(a.tot_npoints), zenang(a.zenang), theta_0(a.theta_0), phi_0(a.phi_0)
{
}
AcxDataSet & AcxDataSet::operator = (AcxDataSet const & a)
{
v_time_data=a.v_time_data; vv_data=a.vv_data; vv_err=a.vv_err;
v_min_data=a.v_min_data; v_max_data=a.v_max_data;
vv_cxdata=a.vv_cxdata; v_min_cxdata=a.v_min_cxdata; v_max_cxdata=a.v_max_cxdata;
tot_npoints=a.tot_npoints; zenang=a.zenang; theta_0=a.theta_0; phi_0=a.phi_0;
return (*this);
}
size_t AcxDataSet::ReadData(TrkInputDataSet & tkds)
{
cout << "---- AcxDataSet::AcxDataSet() reading 4 PAON4 auto-correlation & 6 Cross-cor signals/DataTables for"
<<tkds.NbTrk()<<" tracks ..."<<endl;
if (tkds.NbTrk() != v_time_data.size()) {
v_time_data.resize(tkds.NbTrk());
vv_data.resize(tkds.NbTrk());
vv_err.resize(tkds.NbTrk());
v_min_data.resize(tkds.NbTrk());
v_max_data.resize(tkds.NbTrk());
vv_cxdata.resize(tkds.NbTrk());
vv_cxerr.resize(tkds.NbTrk());
v_min_cxdata.resize(tkds.NbTrk());
v_max_cxdata.resize(tkds.NbTrk());
}
v_freqs=tkds.v_freqs;
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();
tot_npoints = 0; // total number of points for fit
const char * acname[4]={"V11","V22","V33","V44"};
const char * cxname[6]={"V12","V13","V14","V23","V24","V34"};
for(size_t j=0; j<tkds.dataflnm.size(); j++) {
string flnm = tkds.input_base_path+tkds.dataflnm[j]+".ppf";
cout << "1."<<j+1<<" Extracting data from data file DataTable: " << flnm<<endl
<< " ... For time interval (Trk"<<j+1<<") "<<tkds.tstart[j]<<" < t < "<<tkds.tend[j]<<endl;
DataTable dt_data;
PInPersist pin(flnm);
pin >> dt_data;
dt_data.SetShowMinMaxFlag(true);
size_t ktime = dt_data.IndexNom("timesec");
vector<double> vtm;
dt_data.GetColumn(ktime, vtm);
vector< vector<double> > v_vac(NB_ANTENNES);
for(size_t ii=0; ii<NB_ANTENNES; ii++) { // 4 auto-correlations
size_t kac = dt_data.IndexNom(acname[ii]);
dt_data.GetColumn(kac, v_vac[ii]);
vector<double> vtmp, vetmp;
vv_data[j].push_back(vtmp);
vv_err[j].push_back(vetmp);
v_min_data[j].push_back(9.e19);
v_max_data[j].push_back(-9.e19);
}
vector< vector <complex<double> > > v_vcx(NB_CXCORS);
for(size_t ii=0; ii<NB_CXCORS; ii++) { // 6 cross-correlations
size_t kac = dt_data.IndexNom(cxname[ii]);
dt_data.GetColumn(kac, v_vcx[ii]);
vector< complex<double> > vtmp;
vector<double> vetmp;
vv_cxdata[j].push_back(vtmp);
vv_cxerr[j].push_back(vetmp);
v_min_cxdata[j].push_back(9.e19);
v_max_cxdata[j].push_back(-9.e19);
}
vector< vector<double> > & v_data = vv_data[j];
vector< vector<double> > & v_err = vv_err[j];
vector< vector< complex<double> > > & v_cxdata = vv_cxdata[j];
vector< vector<double> > & v_cxerr = vv_cxerr[j];
for(size_t k=0; k<vtm.size(); k++) {
if ((vtm[k]<tkds.tstart[j])||(vtm[k]>tkds.tend[j])) continue;
v_time_data[j].push_back(vtm[k]);
for(size_t ii=0; ii<NB_ANTENNES; ii++) {
vector<double> & vac = v_vac[ii];
v_data[ii].push_back(vac[k]);
v_err[ii].push_back(0.1*sqrt(fabs(vac[k]))); // calcul d'erreur, a affiner
if (vac[k]<v_min_data[j][ii]) v_min_data[j][ii]=vac[k];
if (vac[k]>v_max_data[j][ii]) v_max_data[j][ii]=vac[k];
}
for(size_t ii=0; ii<NB_CXCORS; ii++) { // 6 cross-correlations
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*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;
}
}
tot_npoints += v_time_data[j].size(); // total number of points for fit
cout << " ... Done for " << j+1 << " data size="<<v_time_data[j].size()<<endl;
cout << " Data-AutoCor Min,Max[A1...A4]=";
for(size_t ii=0; ii<NB_ANTENNES; ii++)
cout<<setw(10)<<v_min_data[j][ii]<<","<<setw(10)<<v_max_data[j][ii]<<" ; "; cout << endl;
cout << " Data-CxCorr (abs) Min,Max[Cx1...Cx6]=";
for(size_t ii=0; ii<NB_ANTENNES; ii++)
cout<<setw(10)<<v_min_cxdata[j][ii]<<","<<setw(10)<<v_max_cxdata[j][ii]<<" ; "; cout << endl;
}
return tot_npoints;
}
//------------------------ TrackSet -------------------------------------
TrackSet::TrackSet(TrackSet const & a)
: v_time_sat(a.v_time_sat), v_sat_elev(a.v_sat_elev), v_sat_azim(a.v_sat_azim),
v_interp_elev(a.v_interp_elev), v_interp_sazim(a.v_interp_sazim)
{
}
TrackSet & TrackSet::operator = (TrackSet const & a)
{
v_time_sat=a.v_time_sat; v_sat_elev=a.v_sat_elev; v_sat_azim=a.v_sat_azim;
v_interp_elev=a.v_interp_elev; v_interp_sazim=a.v_interp_sazim;
return *this;
}
TrackSet::TrackSet(TrkInputDataSet & tkds)
{
ReadData(tkds);
}
size_t TrackSet::ReadData(TrkInputDataSet & tkds)
{
cout << "---- TrackSet::ReadData() ; reading source (satellites, ..) for "
<<tkds.NbTrk()<<" tracks ..."<<endl;
if (tkds.NbTrk() != v_time_sat.size()) {
v_time_sat.resize(tkds.NbTrk());
v_sat_elev.resize(tkds.NbTrk());
v_sat_azim.resize(tkds.NbTrk());
v_interp_elev.resize(tkds.NbTrk());
v_interp_sazim.resize(tkds.NbTrk());
}
for(size_t j=0; j<tkds.NbTrk(); j++) {
string flnm = tkds.input_base_path+tkds.trkflnm[j]+".ppf";
cout << "2."<<j+1<<" Extracting data from satellite track DataTables: " << tkds.trkflnm[j] << endl;
DataTable dt_sat;
PInPersist pin(flnm);
pin >> dt_sat;
dt_sat.SetShowMinMaxFlag(true);
size_t ktime = dt_sat.IndexNom("timesec");
dt_sat.GetColumn(ktime, v_time_sat[j]);
size_t kelev = dt_sat.IndexNom("elevation");
dt_sat.GetColumn(kelev, v_sat_elev[j]);
size_t kazim = dt_sat.IndexNom("azimuth");
dt_sat.GetColumn(kazim, v_sat_azim[j]);
cout << "2."<<j+1<<" Done for satellite from file "<<tkds.trkflnm[j]<<" -> size()="<<v_time_sat[j].size()<<endl;
v_interp_elev[j].DefinePoints(v_time_sat[j], v_sat_elev[j]);
double last_azim=v_sat_azim[j][0];
// vector<double> cazim(v_sat_azim[j].size());
// azimuth values, shifted possibly +360 +720 deg ... to avoid jumping from 360 deg to 0 deg
vector<double> shifted_azim(v_sat_azim[j].size());
double azim_offset=0.;
for(size_t k=0; k<v_sat_azim[j].size(); k++) {
double azim=v_sat_azim[j][k];
if ((k>0)&&(azim<last_azim)) {
if ((last_azim>270)&&(azim<90.)) {
azim_offset += 360.;
if (_prtlevel_>0)
cout << " read_srctracks k="<<k<<" last_azim="<<last_azim<<" azimuth= "<<azim<<" Offset->"<<azim_offset<<endl;
}
}
last_azim = azim;
shifted_azim[k]=azim+azim_offset;
/*
double phisrcdeg=90.-v_sat_azim[j][k];
if (phisrcdeg<0.) phisrcdeg+=360.;
double phisrc=Angle(phisrcdeg,Angle::Degree).ToRadian();
cazim[k]=cos(phisrc);
*/
}
v_interp_sazim[j].DefinePoints(v_time_sat[j], shifted_azim);
// v_interp_cphi[j].DefinePoints(v_time_sat[j], cazim);
cout<<" DONE Creation SLinInterp1D for elevation / azimuth ..."<<endl;
cout << v_interp_elev[j];
cout << v_interp_sazim[j];
}
return 0;
}
//------------------------ TrackSet -------------------------------------
ACxSetFitter::ACxSetFitter(AcxDataSet & data, TrackSet & tks)
: fggaussbeam_(true), D_dish(5.), acxd_(data), tks_(tks),
v_RcFit_ac(tks.getNbAutoCor()), v_xi2red_ac(tks.getNbAutoCor()),
v_Ddish(tks.getNbAutoCor()), v_thetaant(tks.getNbAutoCor()),
v_phiant(tks.getNbAutoCor()), v_A(tks.getNbAutoCor()), v_B(tks.getNbAutoCor()),
v_err_Ddish(tks.getNbAutoCor()), v_err_thetaant(tks.getNbAutoCor()),
v_err_phiant(tks.getNbAutoCor()), v_err_A(tks.getNbAutoCor()), v_err_B(tks.getNbAutoCor()),
v_RcFit_cx(tks.getNbCrossCor()), v_xi2red_cx(tks.getNbCrossCor()),
v_phase(tks.getNbCrossCor()), v_Acx(tks.getNbCrossCor()),
v_err_phase(tks.getNbCrossCor()), v_err_Acx(tks.getNbCrossCor())
{
}
int ACxSetFitter::doACfit(string outfilename)
{
cout << "======================================================================================"<<endl;
cout << "---- ACxSetFitter::doACfit() ; Performing antenna pointing fit ..."<<endl;
ofstream ofr(outfilename.c_str());
ofr << "#### Pointing/dish diameter fit on autocorrelation (ACxSetFitter::doACfit() "<<endl
<< "## NumAntenna Xi2red Deff err_Deff Elevation err_Elev Azimuth err_Ezim A0 err_A0 B0 err_B0 A1 err_A1 B1 err_B1 ..."<<endl;
size_t NB_ANTENNES = acxd_.getNbAutoCor();
size_t NTRK = acxd_.NbTrk();
for(size_t ii=0; ii<NB_ANTENNES; ii++) {
v_A[ii].resize(NTRK); v_B[ii].resize(NTRK);
v_err_A[ii].resize(NTRK); v_err_B[ii].resize(NTRK);
}
int tot_npoints_fit = 0;
for(size_t j=0; j<NTRK; j++) tot_npoints_fit += acxd_.v_time_data[j].size();
for(size_t ii=0; ii<NB_ANTENNES; ii++) {
cout << "-------- doACfit() 1."<<ii+1<<" Creating General Fit for AutoCor Antenna= " << ii+1 << endl;
GeneralFitData gdata(1, tot_npoints_fit);
for(size_t j=0; j<NTRK; j++) {
vector< vector<double> > & v_data = acxd_.vv_data[j];
vector< vector<double> > & v_err = acxd_.vv_err[j];
for(size_t k=0; k<acxd_.v_time_data[j].size(); k++) {
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,
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
mFit.SetMaxStep(1000);
// SetParam(int n,double value, double step,double min=1., double max=-1.);
mFit.SetParam(0,"D_dish",D_dish,0.1,D_dish*0.8,D_dish*1.2);
// mFit.SetFix(0, D_dish);
double thetaAntenne=0., phiAntenne=0.;
if (fabs(acxd_.zenang)>1.e-6) {
if (acxd_.zenang<0) {
thetaAntenne=Angle(-acxd_.zenang,Angle::Degree).ToRadian();
phiAntenne=Angle(270.,Angle::Degree).ToRadian();
}
else {
thetaAntenne=Angle(acxd_.zenang,Angle::Degree).ToRadian();
phiAntenne=Angle(90.,Angle::Degree).ToRadian();
}
}
mFit.SetParam(1,"ThetaAntenne",thetaAntenne,M_PI/1440,0.,thetaAntenne+M_PI/36.);
mFit.SetParam(2,"PhiAntenne",phiAntenne,M_PI/180.,0.,2.*M_PI);
// mFit.SetFix(1, thetaAntenne);
// mFit.SetFix(2, phiAntenne);
for(size_t j=0; j<NTRK; j++) {
char pname[32];
sprintf(pname,"A%d",(int)(j+1));
double A = acxd_.v_max_data[j][ii];
mFit.SetParam(2*j+3,pname,A,A/10.,A/20,A*5);
sprintf(pname,"B%d",(int)(j+1));
double B = acxd_.v_min_data[j][ii];
mFit.SetParam(2*j+4,pname,B,B/10.,B/20,B*5);
mFit.SetFix(2*j+4, B);
}
//DBG mFit.PrintFit();
// 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();
v_RcFit_ac[ii]=rcfit; v_xi2red_ac[ii]=-9999.;
if(rcfit>0) {
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()<<" ";
v_xi2red_ac[ii]=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;
ofr <<setw(8)<<Dfit<<" "<<setw(8)<<err_Dfit<<" ";
v_Ddish[ii]=Dfit;
v_err_Ddish[ii]=err_Dfit;
double thetaant=mFit.GetParm(1); double err_thetaant=mFit.GetParmErr(1);
v_thetaant[ii]=thetaant;
double elevdeg=90.-Angle(thetaant).ToDegree();
double err_elevdeg=Angle(err_thetaant).ToDegree();
cout <<setw(16)<<"ThetaAntenne= "<<setw(12)<<Angle(thetaant).ToDegree()<< " +/- "
<<setw(12)<<Angle(err_thetaant).ToDegree()<<" (elevation="
<<setw(8)<<elevdeg<<" +/- "<<setw(8)<<err_elevdeg<<") deg."<<endl;
ofr <<setw(8)<<elevdeg<<" "<<setw(8)<<err_elevdeg<<" ";
double phiant=mFit.GetParm(2); double err_phiant=mFit.GetParmErr(2);
double azimdeg=90.-Angle(phiant).ToDegree();
if (azimdeg<0.) azimdeg += 360.;
double err_azimdeg=Angle(err_phiant).ToDegree();
v_phiant[ii]=phiant;
cout <<setw(16)<<"PhiAntenne= "<<setw(12)<<Angle(phiant).ToDegree()<< " +/- "
<<setw(12)<<Angle(err_phiant).ToDegree()<<" (azimuth ="
<<setw(8)<<azimdeg<<" +/- "<<setw(8)<<err_azimdeg<<" ) deg."<<endl;
ofr <<setw(8)<<azimdeg<<" "<<setw(8)<<err_azimdeg<<" ";
for(size_t j=0; j<NTRK; j++) {
double A=mFit.GetParm(3+2*j); double err_A=mFit.GetParmErr(3+2*j);
double B=mFit.GetParm(4+2*j); double err_B=mFit.GetParmErr(4+2*j);
cout << " Trk/Sat["<<j<<"] -> A= "<<A<<" +/- "<<err_A<<" B= "<<B<<" +/- "<<err_B<<endl;
v_A[ii][j]=A; v_B[ii][j]=B;
v_err_A[ii][j]=err_A; v_err_B[ii][j]=err_B;
ofr <<setw(8)<<A<<" "<<setw(8)<<err_A<<" "<<setw(8)<<B<<" "<<setw(8)<<err_B<<" ";
}
ofr << endl;
}
else {
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);
}
}
return 0;
}
int ACxSetFitter::saveExpectedAC(string outcheckfilename)
{
if (outcheckfilename.length()<1) return 1;
cout << "-----ACxSetFitter::saveExpectedAC() : computing expected signal for fitted params , will be saved to file "
<<outcheckfilename<<endl;
POutPersist pos(outcheckfilename);
size_t NB_ANTENNES = acxd_.getNbAutoCor();
size_t NTRK = acxd_.NbTrk();
for(size_t ii=0; ii<NB_ANTENNES; ii++) {
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,
tks_.v_interp_elev, tks_.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;
}
}
return 0;
}
int ACxSetFitter::doCxfit(string outfilenamecx, string outcheckfilenamecx, bool useAac)
{
size_t NB_ANTENNES=acxd_.getNbAutoCor(); // nombre d'antennes
size_t NB_CXCORS=acxd_.getNbCrossCor();
size_t NTRK = acxd_.NbTrk();
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 ;
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;
int tot_npoints_fit = 0;
for(size_t j=0; j<NTRK; j++) tot_npoints_fit += 2*(acxd_.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 << "--------- 1."<<ii+1<<" doCxfit() Doing fit for CrossCor= " << ii << " FxF= "
<< Anum1[ii]+1<<"x"<<Anum2[ii]+1<<" Baseline="<<baseline<<endl;
GeneralFitData gdata(1, tot_npoints_fit);
for(size_t j=0; j<NTRK; j++) {
vector< vector< complex<double> > > & v_cxdata = acxd_.vv_cxdata[j];
vector< vector<double> > & v_cxerr = acxd_.vv_cxerr[j];
for(size_t k=0; k<acxd_.v_time_data[j].size(); k++) {
gdata.AddData1(acxd_.v_time_data[j][k],v_cxdata[ii][k].real(),v_cxerr[ii][k]); // Fill x, y and error on y
gdata.AddData1(acxd_.v_time_data[j][k],v_cxdata[ii][k].imag(),v_cxerr[ii][k]); // Fill x, y and error on y
}
}
double clight = PhysQty::c().SIValue();
double lambda = clight/(acxd_.v_freqs[0]*1.e6);
ACBeam acb1(v_Ddish[Anum1[ii]], v_thetaant[Anum1[ii]], v_phiant[Anum1[ii]], lambda);
acb1.setGaussianLobe(fggaussbeam_);
ACBeam acb2(v_Ddish[Anum2[ii]], v_thetaant[Anum2[ii]], v_phiant[Anum2[ii]], lambda);
acb2.setGaussianLobe(fggaussbeam_);
CxBeam cxbeam(acb1, acb2, baseline);
MyCxGenXi2 gxi2( acxd_.v_time_data, acxd_.vv_cxdata, acxd_.vv_cxerr,
tks_.v_interp_elev, tks_.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",0.,M_PI/360.,0.,2.2*M_PI);
char oname[32];
vector<double> v_amp(NTRK);
for(size_t j=0; j<NTRK; j++) {
double A=1.; // v_max_cxdata[j][ii];
TVector< complex<double> > signal = gxi2.getDataSignal(j);
Vector asig = SOPHYA::abs(signal);
double mins, maxs;
asig.MinMax(mins, maxs);
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;
/*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.;
double bestxi2 = 9.e19;
double bestphase=0.;
int bestnpts,npts;
int bestafact;
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++) {
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();
double xi2 = gxi2.getXi2(fparm, npts);
if (xi2 < bestxi2) {
bestxi2 = xi2; bestphase=fparm[0]; bestnpts=npts; bestafact=afact[ia];
}
}
}
mFit.SetParam(0,"Phase",bestphase,M_PI/720.,-0.5*M_PI,2.5*M_PI);
cout << "2."<<ii+1<<" Scan param bestxi2_red="<<bestxi2/(double)(tot_npoints_fit-(1+NTRK))<<" 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] << " , ";
double Aac=sqrt(v_A[j][Anum1[ii]] * v_A[j][Anum2[ii]]);
v_Acx[ii][j]=(useAac?Aac:v_amp[j]);
}
cout << endl;
for(size_t j=0; j<NTRK; j++) {
char pname[32];
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));
mFit.SetParam(3+3*j,pname,0.,A/25.,-A/5,A/5.);
mFit.SetFix(2+3*j,0.);
mFit.SetFix(3+3*j,0.);
}
//DBG mFit.PrintFit();
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();
v_RcFit_cx[ii]=rcfit; v_xi2red_cx[ii]=-99999.;
if (_prtlevel_>1) mFit.PrintFit();
if(rcfit>0) {
v_xi2red_cx[ii]=mFit.GetChi2Red();
// 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);
if (phase<0.) phase += 2.*M_PI;
if (phase>2.*M_PI) phase -= 2.*M_PI;
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;
v_err_phase[ii]=err_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["<<j<<"] A= "<<A<<" +/- "<<err_A<<" (A/Ai="<<A/Ai<<")"<<endl;
v_Acx[ii][j]=A; // v_B[ii][j]=B;
v_err_Acx[ii][j]=err_A;
ofr <<setw(8)<<A<<" "<<setw(8)<<err_A<<" ";
}
ofr << endl;
}
else {
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;
}
......@@ -3,13 +3,23 @@
using satellites and celestial sources tracks
R. Ansari, Fevrier 2019 */
#ifndef TRKFIT_H_SEEN
#define TRKFIT_H_SEEN
#include <string>
#include <vector>
#include <iostream>
#include <complex>
#include "slininterp.h"
#include "skymap.h"
#include "histats.h"