Skip to content
GitLab
Projects
Groups
Snippets
Help
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in / Register
Toggle navigation
Open sidebar
BAORadio
AnaPAON4
Commits
25b064ee
Commit
25b064ee
authored
Feb 23, 2019
by
Reza ANSARI
Browse files
Implementation partielle de fit multi-frequence, avec phase=phi0+a*freq, Reza 23/02/2019
parent
1d4158bd
Changes
6
Hide whitespace changes
Inline
Side-by-side
Showing
6 changed files
with
130 additions
and
58 deletions
+130
-58
acbeam.h
acbeam.h
+16
-0
cxbeam.h
cxbeam.h
+21
-0
gcxfit.h
gcxfit.h
+30
-17
makefile
makefile
+3
-12
trkfit.cc
trkfit.cc
+52
-26
trkfit.h
trkfit.h
+8
-3
No files found.
acbeam.h
View file @
25b064ee
...
...
@@ -7,6 +7,7 @@
#include "skymap.h"
#include "unitvector.h"
#include "sunitpcst.h"
using
namespace
std
;
using
namespace
SOPHYA
;
...
...
@@ -90,11 +91,26 @@ public:
{
return
D_
;
}
// define / changes the wavelength for the beam (in meter)
void
setLambda
(
double
lambda
)
{
lambda_
=
lambda
;
DoL_
=
D_
/
lambda
;
return
;
}
// return the associated beam wavelength
inline
double
getLambda
()
const
{
return
lambda_
;
}
// define / changes the wavelength for the beam, by specifying the frequency in MHz
void
setFreq
(
double
freq
)
{
double
clight
=
PhysQty
::
c
().
SIValue
();
double
lambda
=
clight
/
(
freq
*
1.e6
);
return
setLambda
(
lambda
);
}
// Compute the beam normalization factor such as the beam integral over the sphere is equal 1 (unity)
void
AutoNormalizeBeam
(
int
m
=
128
,
int
prtlev
=
0
)
...
...
cxbeam.h
View file @
25b064ee
...
...
@@ -46,6 +46,27 @@ public:
baseline_
=
base_baseline_
;
baseline_
.
Add
(
shift_baseline
);
}
// define / changes the wavelength for the beam (in meter)
void
setLambda
(
double
lambda
)
{
a1_
.
setLambda
(
lambda
);
a2_
.
setLambda
(
lambda
);
facphase_
=
2.
*
M_PI
/
lambda
;
return
;
}
// return the associated beam wavelength
inline
double
getLambda
()
const
{
return
a1_
.
getLambda
();
}
// define / changes the wavelength for the beam, by specifying the frequency in MHz
void
setFreq
(
double
freq
)
{
double
clight
=
PhysQty
::
c
().
SIValue
();
double
lambda
=
clight
/
(
freq
*
1.e6
);
return
setLambda
(
lambda
);
}
//--- return the Beam value for a given direction
inline
complex
<
double
>
Value
(
LongitudeLatitude
const
&
ll
)
{
...
...
gcxfit.h
View file @
25b064ee
...
...
@@ -14,12 +14,12 @@
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_
,
CxBeam
&
beam
_
,
size_t
ncx
=
0
)
:
v_time_data
(
v_time_data_
),
vv_cxdata
(
vv_cxdata_
),
vv_cxerr
(
vv_cxerr_
),
//
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_
),
ncx_
(
ncx
),
beam
(
beam
_
)
ncx_
(
ncx
),
beam
0_
(
beam
),
phlinfac_
(
1.
/
250.
)
{
//DBG cout << " MyCxSignal() CxCorr=ncx "<<ncx_<< endl;
}
...
...
@@ -49,12 +49,20 @@ public:
return
rvs
;
}
virtual
int
getExpectedSignal
(
size_t
j
,
vector
<
complex
<
double
>
>
&
sig
,
double
phase
,
double
A
,
complex
<
double
>
B
=
complex
<
double
>
(
0.
,
0.
))
// Linearly varying phase with frequency
inline
double
getPhase4Freq
(
double
phi0
,
double
a_phi
,
double
freq
)
{
return
(
phi0
+
a_phi
*
(
freq
-
1250.
)
*
phlinfac_
);
}
virtual
int
getExpectedSignal
(
size_t
j
,
vector
<
complex
<
double
>
>
&
sig
,
double
phase
,
double
A
,
complex
<
double
>
B
=
complex
<
double
>
(
0.
,
0.
))
{
//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
];
CxBeam
beam
=
beam0_
;
beam
.
setFreq
(
v_freqs
[
j
]);
// vector< complex<double> > & vcxdata = vv_cxdata[j][nac_];
// vector<double> & vcxerr = vv_cxerr[j][nac_];
for
(
size_t
k
=
0
;
k
<
vtime
.
size
();
k
++
)
{
...
...
@@ -76,9 +84,10 @@ public:
virtual
int
getExpectedSignal
(
size_t
j
,
vector
<
complex
<
double
>
>
&
sig
,
const
double
*
parm
)
{
double
phase
=
parm
[
0
];
double
A
=
parm
[
1
+
3
*
j
];
complex
<
double
>
B
=
complex
<
double
>
(
parm
[
2
+
3
*
j
],
parm
[
3
+
3
*
j
]);
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
);
}
...
...
@@ -116,30 +125,34 @@ public:
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
;
CxBeam
beam
0_
;
size_t
ncx_
;
double
phlinfac_
;
// 1./FreqBandWidth = 1./250.
};
//-----------------------------------------------------------------
class
MyCxGenXi2
:
public
GeneralXi2
,
public
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
<
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
+
3
*
v_time_data_
.
size
())
,
MyCxSignal
(
v_time_data_
,
vv_cxdata_
,
vv_cxerr_
,
v_interp_elev_
,
v_interp_sazim_
,
beam
,
ncx
)
:
GeneralXi2
(
2
+
3
*
v_time_data_
.
size
())
,
MyCxSignal
(
v_time_data_
,
vv_cxdata_
,
vv_cxerr_
,
v_freqs_
,
v_interp_elev_
,
v_interp_sazim_
,
beam
,
ncx
)
{
}
virtual
double
Value
(
GeneralFitData
&
data
,
double
*
parm
,
int
&
ndataused
)
/*
param[0] = phase
parm[1+3*j] = Aj (Amplitude) for track j
parm[2+3*j] = Bj-real Offset, real part
parm[3+3*j] = Bj-imag Offset, imag part
phase = phi_0 + a_phi*(freq-1250.)/250. freq in MHz
param[0] = phi_0
param[1] = a_phi
parm[2+3*j] = Aj (Amplitude) for track j
parm[3+3*j] = Bj-real Offset, real part
parm[4+3*j] = Bj-imag Offset, imag part
*/
{
return
getXi2
(
parm
,
ndataused
);
...
...
makefile
View file @
25b064ee
...
...
@@ -12,13 +12,13 @@ MYOLISTHERE = $(OBJ)/p4autils.o $(OBJ)/visip4reader.o $(OBJ)/visp4winreader.o
MYOLISTHEREFIT
=
$(OBJ)
/trkfit.o
# Define our target list
all
:
rdvisip4 visi2ntac visi2dtacx visi2tmfreq tstutls p4conv2fits msvis2dt visiavg filt_blind ckconvphases vis2ra p4vdblist rdthermrfi visamm
trkacfit
tfm2dt trkacxfit
all
:
rdvisip4 visi2ntac visi2dtacx visi2tmfreq tstutls p4conv2fits msvis2dt visiavg filt_blind ckconvphases vis2ra p4vdblist rdthermrfi visamm tfm2dt trkacxfit
clean
:
rm
-f
$(EXE)
/rdvisip4
$(EXE)
/visi2ntac
$(EXE)
/visi2dtacx
$(EXE)
/visi2tmfreq
$(EXE)
/p4conv2fits
$(EXE)
/msvis2dt
$(EXE)
/visiavg
\
$(EXE)
/filt_blind
$(EXE)
/ckconvphases
$(EXE)
/vis2ra
$(EXE)
/p4vdblist
$(EXE)
/rdthermrfi
$(EXE)
/visamm
$(EXE)
/t
rkacfi
t
$(EXE)
/t
fm2d
t
$(EXE)
/filt_blind
$(EXE)
/ckconvphases
$(EXE)
/vis2ra
$(EXE)
/p4vdblist
$(EXE)
/rdthermrfi
$(EXE)
/visamm
$(EXE)
/t
fm2d
t
$(EXE)
/t
rkacxfi
t
rm
-f
$(OBJ)
/rdvisip4.o
$(OBJ)
/visi2ntac.o
$(OBJ)
/visi2dtacx.o
$(OBJ)
/visi2tmfreq.o
$(OBJ)
/p4conv2fits.o
$(OBJ)
/msvis2dt.o
$(OBJ)
/visiavg.o
\
$(OBJ)
/filt_blind.o
$(OBJ)
/ckconvphases.o
$(OBJ)
/vis2ra.o
$(OBJ)
/p4vdblist.o
$(OBJ)
/rdthermrfi.o
$(OBJ)
/visamm.o
$(OBJ)
/t
rkacfi
t.o
$(OBJ)
/t
fm2d
t.o
$(OBJ)
/filt_blind.o
$(OBJ)
/ckconvphases.o
$(OBJ)
/vis2ra.o
$(OBJ)
/p4vdblist.o
$(OBJ)
/rdthermrfi.o
$(OBJ)
/visamm.o
$(OBJ)
/t
fm2d
t.o
$(OBJ)
/t
rkacxfi
t.o
rm
-f
$(MYOLISTHERE)
depend
:
...
...
@@ -86,15 +86,6 @@ $(OBJ)/visamm.o : visamm.cc $(MYINCLISTHERE)
######
## programme pour fitter l'angle de visee des antennes a partir des autocorrelations + satellites + sources
trkacfit
:
$(EXE)/trkacfit
echo
'---trkacfit made'
$(EXE)/trkacfit
:
$(OBJ)/trkacfit.o $(MYOLISTHERE) $(MYOLISTHEREFIT)
$(CXXLINK)
-o
$(EXE)
/trkacfit
$(OBJ)
/trkacfit.o
$(MYOLISTHERE)
$(MYOLISTHEREFIT)
$(SOPHYAEXTSLBLIST)
$(OBJ)/trkacfit.o
:
trkacfit.cc $(MYINCLISTHERE) $(MYINCLISTHEREFIT)
$(CXXCOMPILE)
-o
$(OBJ)
/trkacfit.o trkacfit.cc
trkacxfit
:
$(EXE)/trkacxfit
echo
'---trkacxfit made'
...
...
trkfit.cc
View file @
25b064ee
...
...
@@ -153,7 +153,7 @@ AcxDataSet::AcxDataSet(AcxDataSet const & a)
tot_npoints
(
a
.
tot_npoints
),
v_freqs
(
a
.
v_freqs
),
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_Acx
(
a
.
v_Acx
),
v_Bcx
(
a
.
v_Bcx
)
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
)
{
}
...
...
@@ -166,7 +166,7 @@ AcxDataSet & AcxDataSet::operator = (AcxDataSet const & a)
tot_npoints
=
a
.
tot_npoints
;
v_freqs
=
a
.
v_freqs
;
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_Acx
=
a
.
v_Acx
;
v_Bcx
=
a
.
v_Bcx
;
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
;
return
(
*
this
);
}
...
...
@@ -370,8 +370,10 @@ ACxSetFitter::ACxSetFitter(AcxDataSet & data, TrackSet & tks)
v_err_phiant
(
tks
.
getNbAutoCor
()),
v_err_A
(
tks
.
getNbAutoCor
()),
v_err_B
(
tks
.
getNbAutoCor
()),
v_acbeams
(
tks
.
getNbAutoCor
()),
v_RcFit_cx
(
tks
.
getNbCrossCor
()),
v_xi2red_cx
(
tks
.
getNbCrossCor
()),
v_phase
(
tks
.
getNbCrossCor
()),
v_Acx
(
tks
.
getNbCrossCor
()),
v_Bcx
(
tks
.
getNbCrossCor
()),
v_err_phase
(
tks
.
getNbCrossCor
()),
v_err_Acx
(
tks
.
getNbCrossCor
()),
v_err_Bcx
(
tks
.
getNbCrossCor
()),
v_phase
(
tks
.
getNbCrossCor
()),
v_phi_0
(
tks
.
getNbCrossCor
()),
v_a_phi
(
tks
.
getNbCrossCor
()),
v_Acx
(
tks
.
getNbCrossCor
()),
v_Bcx
(
tks
.
getNbCrossCor
()),
v_err_phi_0
(
tks
.
getNbCrossCor
()),
v_err_a_phi
(
tks
.
getNbCrossCor
()),
v_err_Acx
(
tks
.
getNbCrossCor
()),
v_err_Bcx
(
tks
.
getNbCrossCor
()),
v_cxbeams
(
tks
.
getNbCrossCor
())
{
if
(
data
.
NbTrk
()
!=
tks
.
NbTrk
())
...
...
@@ -547,7 +549,7 @@ int ACxSetFitter::saveExpectedAC(string outcheckfilename)
}
int
ACxSetFitter
::
doCxfit
(
string
outfilenamecx
,
bool
useAac
)
int
ACxSetFitter
::
doCxfit
(
string
outfilenamecx
,
bool
useAac
,
bool
fgphi0only
)
{
size_t
NB_ANTENNES
=
acxd_
.
getNbAutoCor
();
// nombre d'antennes
size_t
NB_CXCORS
=
acxd_
.
getNbCrossCor
();
...
...
@@ -558,7 +560,7 @@ int ACxSetFitter::doCxfit(string outfilenamecx, bool useAac)
if
(
useAac
)
cout
<<
" ... Using Amplitude from auto-correlations fit for initial fit parameter value..."
<<
endl
;
ofstream
ofr
(
outfilenamecx
.
c_str
());
ofr
<<
"#### cross-cor phase fit (ACxSetFitter::doCxfit() ) "
<<
endl
<<
"## NumCxCor RcFit Xi2red
Phase
err_Ph
ase
(deg)
A0 err_A0 A1 err_A1 ..."
<<
endl
;
<<
"## NumCxCor RcFit Xi2red
Phi0 err_Phi0 a_Phi
err_
a_
Ph
i
(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
());
cout
<<
" Total number of data points for fit="
<<
tot_npoints_fit
<<
endl
;
...
...
@@ -595,13 +597,19 @@ int ACxSetFitter::doCxfit(string outfilenamecx, bool useAac)
CxBeam
cxbeam
(
acb1
,
acb2
,
baseline
);
v_cxbeams
[
ii
]
=
cxbeam
;
MyCxGenXi2
gxi2
(
acxd_
.
v_time_data
,
acxd_
.
vv_cxdata
,
acxd_
.
vv_cxerr
,
MyCxGenXi2
gxi2
(
acxd_
.
v_time_data
,
acxd_
.
vv_cxdata
,
acxd_
.
vv_cxerr
,
acxd_
.
v_freqs
,
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
(
1
000
);
mFit
.
SetMaxStep
(
3
000
);
// 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
);
mFit
.
SetParam
(
0
,
"Phi_0"
,
0.
,
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
;
mFit
.
SetFix
(
1
,
0.
);
}
else
cout
<<
" ACxSetFitter::doCxfit() Fitting Phase(freq) = Phi0 + a_Phi * (freq-1250.)/250. "
<<
endl
;
char
oname
[
32
];
vector
<
double
>
v_amp
(
NTRK
);
...
...
@@ -629,21 +637,22 @@ int ACxSetFitter::doCxfit(string outfilenamecx, bool useAac)
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.
;
fparm
[
2
+
3
*
j
]
=
(
useAac
?
Aac
:
v_amp
[
j
]);
fparm
[
2
+
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
();
fparm
[
1
]
=
0.
;
double
xi2
=
gxi2
.
getXi2
(
fparm
,
npts
);
if
(
xi2
<
bestxi2
)
{
bestxi2
=
xi2
;
bestphase
=
fparm
[
0
];
bestnpts
=
npts
;
bestafact
=
afact
[
ia
];
}
}
}
mFit
.
SetParam
(
0
,
"Ph
ase
"
,
bestphase
,
M_PI
/
720.
,
-
0.5
*
M_PI
,
2.5
*
M_PI
);
mFit
.
SetParam
(
0
,
"Ph
i_0
"
,
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_ph
ase
[
ii
]
=
bestphase
;
v_ph
i_0
[
ii
]
=
bestphase
;
for
(
size_t
j
=
0
;
j
<
NTRK
;
j
++
)
{
cout
<<
v_amp
[
j
]
<<
" , "
;
double
Aac
=
sqrt
(
v_A
[
Anum1
[
ii
]][
j
]
*
v_A
[
Anum2
[
ii
]][
j
]);
...
...
@@ -656,13 +665,13 @@ int ACxSetFitter::doCxfit(string outfilenamecx, bool useAac)
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
);
mFit
.
SetParam
(
2
+
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.
);
sprintf
(
pname
,
"Bim%d"
,(
int
)(
j
+
1
));
mFit
.
SetParam
(
4
+
3
*
j
,
pname
,
0.
,
A
/
25.
,
-
A
/
5
,
A
/
5.
);
mFit
.
SetFix
(
3
+
3
*
j
,
0.
);
mFit
.
SetFix
(
4
+
3
*
j
,
0.
);
}
//DBG mFit.PrintFit();
if
(
_prtlevel_
>
1
)
...
...
@@ -683,17 +692,27 @@ int ACxSetFitter::doCxfit(string outfilenamecx, bool useAac)
}
ofr
<<
setw
(
4
)
<<
ii
+
1
<<
" "
<<
setw
(
5
)
<<
rcfit
<<
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
()
<<
" "
;
double
phi0
=
mFit
.
GetParm
(
0
);
double
err_phi0
=
mFit
.
GetParmErr
(
0
);
double
aphi
=
mFit
.
GetParm
(
1
);
double
err_aphi
=
mFit
.
GetParmErr
(
1
);
// on calcule la phase ajustee pour la frequence de reference 1300 MHz
double
phase
=
gxi2
.
getPhase4Freq
(
phi0
,
aphi
,
1300.
);
while
(
phase
<
0.
)
phase
+=
2.
*
M_PI
;
while
(
phase
>
2.
*
M_PI
)
phase
-=
2.
*
M_PI
;
cout
<<
"Phase(@1300MHz)= "
<<
setw
(
10
)
<<
Angle
(
phase
).
ToDegree
()
<<
" phi_0= "
<<
setw
(
10
)
<<
Angle
(
phi0
).
ToDegree
()
<<
" +/- "
<<
setw
(
10
)
<<
Angle
(
err_phi0
).
ToDegree
()
<<
" deg."
<<
" a_phi= "
<<
setw
(
8
)
<<
Angle
(
aphi
).
ToDegree
()
<<
" +/- "
<<
setw
(
10
)
<<
Angle
(
err_aphi
).
ToDegree
()
<<
" deg/250 MHz"
<<
endl
;
ofr
<<
setw
(
8
)
<<
Angle
(
phi0
).
ToDegree
()
<<
" "
<<
setw
(
8
)
<<
Angle
(
err_phi0
).
ToDegree
()
<<
" "
<<
setw
(
8
)
<<
Angle
(
aphi
).
ToDegree
()
<<
" "
<<
setw
(
8
)
<<
Angle
(
err_aphi
).
ToDegree
()
<<
" "
;
v_phase
[
ii
]
=
phase
;
v_err_phase
[
ii
]
=
err_phase
;
v_phi_0
[
ii
]
=
phi0
;
v_err_phi_0
[
ii
]
=
err_phi0
;
v_a_phi
[
ii
]
=
aphi
;
v_err_a_phi
[
ii
]
=
err_aphi
;
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
);
double
A
=
mFit
.
GetParm
(
2
+
3
*
j
);
double
err_A
=
mFit
.
GetParmErr
(
2
+
3
*
j
);
cout
<<
" Trk["
<<
j
<<
"] A= "
<<
A
<<
" +/- "
<<
err_A
<<
" (A/Ai="
<<
A
/
Ai
<<
")"
<<
endl
;
v_Acx
[
ii
][
j
]
=
A
;
v_Bcx
[
ii
][
j
]
=
complex
<
double
>
(
0.
,
0.
);
...
...
@@ -716,6 +735,8 @@ int ACxSetFitter::doCxfit(string outfilenamecx, bool useAac)
fit_cx_done
=
true
;
acxd_
.
v_cxbeams
=
v_cxbeams
;
acxd_
.
v_phase
=
v_phase
;
acxd_
.
v_phi_0
=
v_phi_0
;
acxd_
.
v_a_phi
=
v_a_phi
;
acxd_
.
v_Acx
=
v_Acx
;
acxd_
.
v_Bcx
=
v_Bcx
;
return
0
;
...
...
@@ -733,13 +754,15 @@ int ACxSetFitter::saveExpectedCx(string outcheckfilename)
for
(
size_t
ii
=
0
;
ii
<
NB_CXCORS
;
ii
++
)
{
CxBeam
cxbeam
=
v_cxbeams
[
ii
];
MyCxSignal
cxsig
(
acxd_
.
v_time_data
,
acxd_
.
vv_cxdata
,
acxd_
.
vv_cxerr
,
MyCxSignal
cxsig
(
acxd_
.
v_time_data
,
acxd_
.
vv_cxdata
,
acxd_
.
vv_cxerr
,
acxd_
.
v_freqs
,
tks_
.
v_interp_elev
,
tks_
.
v_interp_sazim
,
cxbeam
,
ii
);
for
(
size_t
j
=
0
;
j
<
NTRK
;
j
++
)
{
TVector
<
complex
<
double
>
>
signal
=
cxsig
.
getDataSignal
(
j
);
sprintf
(
oname
,
"cx_%d_%d"
,(
int
)
ii
+
1
,(
int
)
j
+
1
);
pox
<<
PPFNameTag
(
oname
)
<<
signal
;
TVector
<
complex
<
double
>
>
expsignal
=
cxsig
.
getExpectedSignal
(
j
,
v_phase
[
ii
],
v_Acx
[
ii
][
j
]);
//DBG cout << " *DBG* getPhase4Freq() phi0="<<acxd_.v_phi_0[ii]<<" a_phi="<<acxd_.v_a_phi[ii]<<" freq="<<acxd_.v_freqs[j]<<endl;
double
phase
=
cxsig
.
getPhase4Freq
(
acxd_
.
v_phi_0
[
ii
],
acxd_
.
v_a_phi
[
ii
],
acxd_
.
v_freqs
[
j
]);
TVector
<
complex
<
double
>
>
expsignal
=
cxsig
.
getExpectedSignal
(
j
,
phase
,
v_Acx
[
ii
][
j
]);
sprintf
(
oname
,
"simcx_%d_%d"
,(
int
)
ii
+
1
,(
int
)
j
+
1
);
pox
<<
PPFNameTag
(
oname
)
<<
expsignal
;
if
(
ii
==
0
)
{
...
...
@@ -783,6 +806,7 @@ CxBaselineFitter::~CxBaselineFitter()
void
CxBaselineFitter
::
initFitParams
()
{
//DBG cout << " *DBG* CxBaselineFitter::initFitParams() v_acxd[0].v_phase.size()="<<v_acxd[0].v_phase.size()<<endl;
size_t
NB_ANTENNES
=
v_acxd
[
0
].
getNbAutoCor
();
// nombre d'antennes
size_t
NB_CXCORS
=
v_acxd
[
0
].
getNbCrossCor
();
if
(
NB_ANTENNES
!=
4
)
...
...
@@ -801,6 +825,8 @@ void CxBaselineFitter::initFitParams()
bestfitparam
[
3
*
(
i
+
1
)
+
j
]
=
err_bestfitparam
[
3
*
(
i
+
1
)
+
j
]
=
0.
;
}
}
//DBG cout << " *DBG* DONE **** CxBaselineFitter::initFitParams()"<<endl;
}
int
CxBaselineFitter
::
dofit
(
string
outfilename
,
bool
fgfixbaseline
)
...
...
trkfit.h
View file @
25b064ee
...
...
@@ -86,7 +86,9 @@ public:
vector
<
CxBeam
>
v_acbeams
;
// the four aut-correlations beams after fit
vector
<
CxBeam
>
v_cxbeams
;
// the six cross correlations after fit
vector
<
double
>
v_phase
;
// fitted phases for individual cross-correlations
vector
<
double
>
v_phase
;
// fitted phases for individual cross-correlations at reference frequency
vector
<
double
>
v_phi_0
;
// fitted phases (Phi0) for individual cross-correlations
vector
<
double
>
v_a_phi
;
// fitted phases (a_phi = slope with frequency) for individual cross-correlations
vector
<
vector
<
double
>
>
v_Acx
;
// Amplitudes for the six cross-cors after fit
vector
<
vector
<
complex
<
double
>
>
>
v_Bcx
;
// Offset for the six cross-cors after fit
};
...
...
@@ -128,7 +130,7 @@ class ACxSetFitter {
public:
ACxSetFitter
(
AcxDataSet
&
data
,
TrackSet
&
tks
);
int
doACfit
(
string
outfilename
);
// outfilename : fit results
int
doCxfit
(
string
outfilenamecx
,
bool
useAac
=
true
);
// outfilename : fit results
int
doCxfit
(
string
outfilenamecx
,
bool
useAac
=
true
,
bool
fgphi0only
=
true
);
// outfilename : fit results
int
saveExpectedAC
(
string
outcheckfilename
);
// outfilename PPF (or fits ?) file with the expected AC signals
int
saveExpectedCx
(
string
outcheckfilename
);
// outfilename PPF (or fits ?) file with the expected signals
...
...
@@ -154,9 +156,12 @@ public:
vector
<
int
>
v_RcFit_cx
;
vector
<
double
>
v_xi2red_cx
;
vector
<
double
>
v_phase
;
vector
<
double
>
v_phi_0
;
vector
<
double
>
v_a_phi
;
vector
<
vector
<
double
>
>
v_Acx
;
vector
<
vector
<
complex
<
double
>
>
>
v_Bcx
;
// Offset for the six cross-cors after fit
vector
<
double
>
v_err_phase
;
vector
<
double
>
v_err_phi_0
;
vector
<
double
>
v_err_a_phi
;
vector
<
vector
<
double
>
>
v_err_Acx
;
vector
<
vector
<
complex
<
double
>
>
>
v_err_Bcx
;
vector
<
CxBeam
>
v_cxbeams
;
// the six cross correlations after fit
...
...
Write
Preview
Markdown
is supported
0%
Try again
or
attach a new file
.
Attach a file
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment