Skip to content

GitLab

Commit 9aa276e6 authored by Reza ANSARI's avatar Reza ANSARI
Browse files

debug des fits avec phase(freq), Reza 27/02/2019

parent adc5825d
Hide whitespace changes
Inline Side-by-side
Showing with 27 additions and 15 deletions
gacfit.h
View file @ 9aa276e6
......@@ -114,8 +114,8 @@ public:
double B = 1.;
if (v_idxparam[j]>=0) {
int jj = v_idxparam[j];
double A = parm[3+2*jj];
double B = parm[4+2*jj];
A = parm[3+2*jj];
B = parm[4+2*jj];
}
return getExpectedSignal(j, sig, D_dish, theta_b, phi_b, A, B);
}
......
gcxfit.h
View file @ 9aa276e6
......@@ -81,8 +81,8 @@ public:
complex<double> cxb=beam.Value(thetasrc, Angle(phisrcdeg,Angle::Degree).ToRadian());
// 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;
//DBG cout << " *DBG* "<<k<<" Theta,Phi-src = " << Angle(thetasrc).ToDegree() << " , " << phisrcdeg
// << " elev="<<elev<<" --> cxb= "<<cxb<< " sig="<<sig[k]<<" A="<<A<<" phase="<<phase<<endl;
}
return (int)sig.size();
......@@ -93,8 +93,9 @@ public:
double phase = getPhase4Freq(parm[0], parm[1], v_freqs[j]);
double A = parm[2+3*j];
complex<double> B = complex<double>(parm[3+3*j], parm[4+3*j]);
return getExpectedSignal(j, sig, phase, A);
//DBG cout << " getExpectedSignal()*DBG* : parm[0...2]:"<<parm[0]<<" "<<parm[1]<<" "<<parm[2]<<" -> Ph="
//DBG << phase<<" A="<<A<<endl;
return getExpectedSignal(j, sig, phase, A, B);
}
virtual TVector< complex<double> > getExpectedSignal(size_t j, double phase, double A, complex<double> B=complex<double>(0.,0.)) const
......
trkfit.cc
View file @ 9aa276e6
......@@ -434,11 +434,15 @@ int ACxSetFitter::doACfit(string outfilename)
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 j=0; j<NTRK; j++) {
if (acxd_.v_noAC[j]) continue;
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++) {
if (acxd_.v_noAC[j]) continue;
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++) {
......@@ -450,9 +454,9 @@ int ACxSetFitter::doACfit(string outfilename)
// GeneralFit mFit(&gxi2);
TkF_Fitter mFit(gxi2);
mFit.SetData(&gdata); // connect data to the fitter , here the data is unused - gxi2 includes its data
mFit.SetMaxStep(1000);
mFit.SetMaxStep(5000);
// 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.SetParam(0,"D_dish",D_dish,0.1,D_dish*0.7,D_dish*1.4);
// mFit.SetFix(0, D_dish);
double thetaAntenne=0., phiAntenne=0.;
......@@ -466,22 +470,25 @@ int ACxSetFitter::doACfit(string outfilename)
phiAntenne=Angle(90.,Angle::Degree).ToRadian();
}
}
mFit.SetParam(1,"ThetaAntenne",thetaAntenne,M_PI/1440,0.,M_PI/3.);
mFit.SetParam(1,"ThetaAntenne",thetaAntenne,M_PI/1440,0.,M_PI/4.); // thetaAntenne+M_PI/30.); //
mFit.SetParam(2,"PhiAntenne",phiAntenne,M_PI/180.,0.,2.*M_PI);
// mFit.SetFix(1, thetaAntenne);
// mFit.SetFix(2, phiAntenne);
size_t jj=0;
for(size_t j=0; j<NTRK; j++) {
double A = acxd_.v_max_data[j][ii];
double B = acxd_.v_min_data[j][ii];
v_A[ii][j]=A; v_err_A[ii][j]=0.;
v_B[ii][j]=B; v_err_B[ii][j]=0.;
if (acxd_.v_noAC[j]) continue;
char pname[32];
sprintf(pname,"A%d",(int)(j+1));
double A = acxd_.v_max_data[j][ii];
mFit.SetParam(2*jj+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*jj+4,pname,B,B/10.,B/20,B*5);
// mFit.SetFix(2*jj+4, B);
mFit.SetFix(2*jj+4, B);
jj++;
}
//DBG mFit.PrintFit();
// cout << "do_p4_trkfit 2."<<ii+1<<" Performing the fit for AutoCor Antenna= " << ii+1 << endl;
......@@ -533,6 +540,7 @@ int ACxSetFitter::doACfit(string outfilename)
B=mFit.GetParm(4+2*jj); err_B=mFit.GetParmErr(4+2*jj);
jj++;
}
v_A[ii][j]=A; v_err_A[ii][j]=err_A; v_B[ii][j]=B; v_err_B[ii][j]=err_B;
cout << " Trk/Sat["<<j<<"] -> A= "<<A<<" +/- "<<err_A<<" B= "<<B<<" +/- "<<err_B<<(acxd_.v_noAC[j]?" NOT Fitted":"")<<endl;
if (acxd_.v_noAC[j])
ofr <<setw(8)<<A<<" "<<setw(8)<<" NOFIT "<<" "<<setw(8)<<B<<" "<<setw(8)<<" NOFIT "<<" ";
......@@ -652,7 +660,7 @@ int ACxSetFitter::doCxfit(string outfilenamecx, bool useAac, bool fgphi0only)
mFit.SetData(&gdata); // connect data to the fitter , here the data is unused - gxi2 includes its data
mFit.SetMaxStep(3000);
// SetParam(int n,double value, double step,double min=1., double max=-1.);
mFit.SetParam(0,"Phi_0",0.,M_PI/360.,0.,2.2*M_PI);
mFit.SetParam(0,"Phi_0",0.5*M_PI,M_PI/360.,0.,2.2*M_PI);
mFit.SetParam(1,"a_phi",0.,0.05,-15.,15.);
if (fgphi0only) {
cout << " ACxSetFitter::doCxfit() Fitting Phi0 Only (frequency independent phase)"<<endl;
......@@ -686,13 +694,16 @@ int ACxSetFitter::doCxfit(string outfilenamecx, bool useAac, bool fgphi0only)
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]);
//DBG cout << " *DBG* j="<<j<<" ia="<<ia<<" vA="<<v_A[Anum1[ii]][j]<<" x "<<v_A[Anum2[ii]][j]
// <<" -> "<<Aac<<endl;
fparm[2+3*j]=(useAac?Aac:v_amp[j]);
fparm[2+3*j]*=afact[ia]; fparm[2+3*j]=fparm[3+3*j]=0.;
fparm[2+3*j]*=afact[ia]; fparm[3+3*j]=fparm[4+3*j]=0.;
}
for(double ph=0.; ph<360.; ph += 1) {
fparm[0]=Angle(ph, Angle::Degree).ToRadian();
fparm[1]=0.;
double xi2 = gxi2.getXi2(fparm, npts);
//DBG cout << " *DBG* ia="<<ia<<" afact="<<afact[ia]<<" ph="<<ph<<" xi2="<<xi2<<endl;
if (xi2 < bestxi2) {
bestxi2 = xi2; bestphase=fparm[0]; bestnpts=npts; bestafact=afact[ia];
}
......
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