Skip to content

GitLab

Commit 673e0ba0 authored by Reza ANSARI's avatar Reza ANSARI
Browse files

Implementation de la logique pour prise en compte flag NOCX (ne pas prendre en... 

Implementation de la logique pour prise en compte flag NOCX (ne pas prendre en compte certaines traces pour le fit des cross-cor, Reza 03/03/2019
parent e998f1eb
Hide whitespace changes
Inline Side-by-side
Showing with 25 additions and 13 deletions
gcxfit.h
View file @ 673e0ba0
......@@ -16,10 +16,10 @@
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<bool> & v_noCx_,
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_),
: v_time_data(v_time_data_), vv_cxdata(vv_cxdata_), vv_cxerr(vv_cxerr_), v_freqs(v_freqs_), v_noCx(v_noCx_),
v_interp_elev(v_interp_elev_), v_interp_sazim(v_interp_sazim_),
ncx_(ncx), beam0_(beam), phlinfac_(1./250.)
{
......@@ -114,6 +114,7 @@ public:
double xi2=0.;
npts=0;
for(size_t j=0; j<v_time_data.size(); j++) {
if (v_noCx[j]) continue; // pour sauter les signaux flagges a ne pas etre utilise en Cross-Corr
vector< complex<double> > sig;
npts += getExpectedSignal(j, sig, parm);
vector< complex<double> > & vcxdata = vv_cxdata[j][ncx_];
......@@ -133,6 +134,7 @@ public:
vector< vector< vector< complex<double> > > > & vv_cxdata;
vector< vector< vector<double> > > & vv_cxerr;
vector<double> & v_freqs;
vector<bool> & v_noCx;
vector< SLinInterp1D > & v_interp_elev;
vector< SLinInterp1D > & v_interp_sazim;
CxBeam beam0_;
......@@ -145,11 +147,11 @@ public:
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<bool> & v_noCx_,
vector< SLinInterp1D > & v_interp_elev_, vector< SLinInterp1D > & v_interp_sazim_,
CxBeam & beam, size_t ncx=0)
: 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)
MyCxSignal(v_time_data_, vv_cxdata_, vv_cxerr_, v_freqs_, v_noCx_, v_interp_elev_, v_interp_sazim_, beam, ncx)
{
}
......@@ -177,10 +179,10 @@ typedef MyCxGenXi2 TkF_CxXi2 ;
class MyCx_Xi2_Minuit : public Trk_FCNBase , public MyCxSignal {
public:
MyCx_Xi2_Minuit(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<bool> & v_noCx_,
vector< SLinInterp1D > & v_interp_elev_, vector< SLinInterp1D > & v_interp_sazim_,
CxBeam & beam, size_t ncx=0)
: MyCxSignal(v_time_data_, vv_cxdata_, vv_cxerr_, v_freqs_, v_interp_elev_, v_interp_sazim_, beam, ncx), errDef_(1.)
: MyCxSignal(v_time_data_, vv_cxdata_, vv_cxerr_, v_freqs_, v_noCx_, v_interp_elev_, v_interp_sazim_, beam, ncx), errDef_(1.)
{
nbparam_= 2+3*v_time_data_.size();
myparam_ = new double[nbparam_];
......
trkfit.cc
View file @ 673e0ba0
......@@ -619,7 +619,7 @@ int ACxSetFitter::doCxfit(string outfilenamecx, bool useAac, bool fg_B0, bool fg
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 Phi0 err_Phi0 a_Phi err_a_Phi (deg) A0 err_A0 A1 err_A1 ..."<<endl;
<< "## NumCxCor RcFit Xi2red Phi0 err_Phi0 a_Phi err_a_Phi (deg) A0 err_A0 B0 errB0 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;
......@@ -656,7 +656,7 @@ int ACxSetFitter::doCxfit(string outfilenamecx, bool useAac, bool fg_B0, bool fg
CxBeam cxbeam(acb1, acb2, baseline);
v_cxbeams[ii]=cxbeam;
TkF_CxXi2 gxi2( acxd_.v_time_data, acxd_.vv_cxdata, acxd_.vv_cxerr, acxd_.v_freqs,
TkF_CxXi2 gxi2( acxd_.v_time_data, acxd_.vv_cxdata, acxd_.vv_cxerr, acxd_.v_freqs, acxd_.v_noCx,
tks_.v_interp_elev, tks_.v_interp_sazim, cxbeam, ii); // MyCxGenXi2
// GeneralFit mFit(&gxi2);
TkF_Fitter mFit(gxi2);
......@@ -733,10 +733,17 @@ int ACxSetFitter::doCxfit(string outfilenamecx, bool useAac, bool fg_B0, bool fg
mFit.SetParam(3+3*j,pname,0.,A/25.,-A/5,A/5.);
sprintf(pname,"Bim%d",(int)(j+1));
mFit.SetParam(4+3*j,pname,0.,A/25.,-A/5,A/5.);
if (fg_B0) {
if (acxd_.v_noCx[j]) {
mFit.SetFix(2+3*j,A);
mFit.SetFix(3+3*j,0.);
mFit.SetFix(4+3*j,0.);
}
else {
if (fg_B0) {
mFit.SetFix(3+3*j,0.);
mFit.SetFix(4+3*j,0.);
}
}
}
//DBG mFit.PrintFit();
if (_prtlevel_>1)
......@@ -778,11 +785,14 @@ int ACxSetFitter::doCxfit(string outfilenamecx, bool useAac, bool fg_B0, bool fg
double Aac=sqrt(v_A[Anum1[ii]][j] * v_A[Anum2[ii]][j]);
double Ai=(useAac?Aac:v_amp[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;
complex<double> B(mFit.GetParm(3+3*j), mFit.GetParm(4+3*j));
complex<double> err_B(mFit.GetParmErr(3+3*j), mFit.GetParmErr(4+3*j));
cout << " Trk["<<j<<"] A= "<<A<<" +/- "<<err_A<<" (A/Ai="<<A/Ai<<")"<<
" B= "<<B<<" +/- "<<err_B<<(acxd_.v_noCx[j]?" NoFIT":" ")<<endl;
v_Acx[ii][j]=A;
v_Bcx[ii][j]=complex<double>(0.,0.);
v_Bcx[ii][j]=B;
v_err_Acx[ii][j]=err_A;
ofr <<setw(8)<<A<<" "<<setw(8)<<err_A<<" ";
ofr <<setw(8)<<A<<" "<<setw(8)<<err_A<<" "<<setw(14)<<B<<" "<<setw(12)<<err_B<<" ";
}
ofr << endl;
}
......@@ -819,7 +829,7 @@ 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, acxd_.v_freqs,
MyCxSignal cxsig( acxd_.v_time_data, acxd_.vv_cxdata, acxd_.vv_cxerr, acxd_.v_freqs, acxd_.v_noCx,
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);
......
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