Commit 48d5e7b3 authored by Reza  ANSARI's avatar Reza ANSARI
Browse files

Implementation classe CxBaselineFitter ds trkfit.h (Ajustement baselines et...

Implementation classe CxBaselineFitter ds trkfit.h (Ajustement baselines et phases) et appel ds trkacxfit.cc , Reza 19/02/2019
parent ba281357
......@@ -88,6 +88,8 @@ static int prtlevel=0; // print level
static string default_input_path;
static vector<string> trksetfiles; // datacard files defining the track/data sets
static bool do_baselinefit = false; // if true , perform baseline fit
//--- End of list of global parameters
......@@ -132,6 +134,13 @@ int main (int narg, char* arg[])
v_acxd.push_back(acxd);
v_trk.push_back(tks);
}
if (do_baselinefit) {
CxBaselineFitter cxbfit(v_acxd, v_trk);
string cxbofile = "cxb6_"+outfilename;
string cxbckfile = "cxb6_"+checkfilename;
cxbfit.dofit(cxbofile, cxbckfile);
}
Timer tm("trkacxfit");
} // End of try bloc
catch (PThrowable & exc) {
......@@ -167,7 +176,7 @@ int decode_args(int narg, char* arg[])
cout << " trkacxfit : fit array geometry and \n"
<< " Usage: trkacxfit [-options] track_Set_1 [track_Set_2] [track_Set_3 ...] \n"
<< " options: -inp def_input_path -out OutFileName -ckf CheckFileName \n"
<< " -D dish_diameter -ngb -prt PrintLevel \n"<<endl;
<< " -doblf -ngb -D dish_diameter -prt PrintLevel \n"<<endl;
if (!fglonghelp) {
cout << " trkacxfit -h to get option description " << endl;
return 2;
......@@ -178,8 +187,9 @@ int decode_args(int narg, char* arg[])
<< " -ckf CheckFileName : Output file (PPF) to save expected signals (def chktrkfit.ppf) \n"
<< " Output file names are constructed from OutFileName and CheckFileName prepending \n"
<< " ac_JJ_ or cx_JJ_ where JJ=1...n is the track_set number \n"
<< " -doblf : Perform baseline determination (fit) - default NO \n"
<< " -ngb : Use non gaussian beam profile (Bessel j1(angle) - default Gaussian beam) \n"
<< " -D dish_diameter : define effective dish diameter (in meter D*eff , def=4.5) \n"
<< " -ngb : Use non gaussian beam profile (Bessel j1(angle)) \n"
<< " -prt PrintLevel: specify print level \n"
<< " track_Set_J (J=1..n) are the auto/cross correlation and track data sets \n"
<< " in datacard format [ @trk (multiple) @zenang cards @inpath] \n"
......@@ -204,6 +214,9 @@ int decode_args(int narg, char* arg[])
else if (fbo=="-ngb") { // Use Non gaussian beam
fggaussbeam=false; arg++; narg--; lastargs.clear();
}
else if (fbo=="-doblf") { // Use Non gaussian beam
do_baselinefit=false; arg++; narg--; lastargs.clear();
}
else if (fbo=="-inp") { // output file name
if (narg<2) { cout << "trkacxfit/decode_args missing/bad argument, -h for help " << endl; return -1; }
default_input_path=arg[2]; arg+=2; narg-=2; lastargs.clear();
......@@ -231,6 +244,7 @@ int decode_args(int narg, char* arg[])
cout << " ------------------- trkacxfit/run parameters:"<<endl;
cout << " Beam: D_dish (initial guess)="<<D_dish<<" GaussianBeam ? "<<(fggaussbeam?"Yes":"No")<<endl;
cout << " OutFileName= "<<outfilename<<" CheckFileName= "<<checkfilename<<endl;
cout << " Perform baseline fit ? " << (do_baselinefit?"Yes":"No")<<endl;
cout << " --- TrackSetFiles: (NbFiles="<<lastargs.size()<<" default directory: "<<default_input_path<<")"<<endl;
trksetfiles = lastargs;
for (size_t i=0; i<trksetfiles.size(); i++) {
......
......@@ -337,7 +337,7 @@ size_t TrackSet::ReadData(TrkInputDataSet & tkds)
}
//------------------------ TrackSet -------------------------------------
//------------------------ ACxSetFitter -------------------------------------
ACxSetFitter::ACxSetFitter(AcxDataSet & data, TrackSet & tks)
: fggaussbeam_(true), D_dish(5.), acxd_(data), tks_(tks), fit_ac_done(false), fit_cx_done(false),
v_RcFit_ac(tks.getNbAutoCor()), v_xi2red_ac(tks.getNbAutoCor()),
......@@ -351,7 +351,10 @@ ACxSetFitter::ACxSetFitter(AcxDataSet & data, TrackSet & tks)
v_err_phase(tks.getNbCrossCor()), v_err_Acx(tks.getNbCrossCor()), v_err_Bcx(tks.getNbCrossCor()),
v_cxbeams(tks.getNbCrossCor())
{
if (data.NbTrk() != tks.NbTrk())
throw ParmError("ACxSetFitter(data, tks) NOT same number of tracks NbTrk() in data and tks");
if (data.NbTrk() < 1)
throw ParmError("ACxSetFitter(data, tks) 0 tracks in data data.NbTrk()<1 ");
}
int ACxSetFitter::doACfit(string outfilename)
......@@ -690,7 +693,90 @@ int ACxSetFitter::doCxfit(string outfilenamecx, string outcheckfilenamecx, bool
if (pox) delete pox;
fit_cx_done=true;
acxd_.v_cxbeams=v_cxbeams;
acxd_.v_phase=v_phase;
acxd_.v_Acx=v_Acx;
acxd_.v_Bcx=v_Bcx;
return 0;
}
//------------------------ CxBaselineFitter -------------------------------------
CxBaselineFitter::CxBaselineFitter(vector<AcxDataSet> & v_data, vector<TrackSet> & v_tks)
: v_acxd(v_data), v_trks(v_tks), tot_ntrks(0), fit_done(false), xi2red(-9.e9)
{
if (v_acxd.size() != v_trks.size())
throw ParmError("CxBaselineFitter::CxBaselineFitter(v_data, v_tks) NOT same size v_data,v_tks ");
if (v_acxd.size() < 1)
throw ParmError("CxBaselineFitter::CxBaselineFitter(v_data, v_tks) v_data.size()<1 ");
v_phases.resize(v_acxd[0].getNbAutoCor()-1);
v_baselineshits.resize(v_acxd[0].getNbAutoCor()-1);
tot_ntrks=0;
for(size_t i=0; i<v_acxd.size(); i++) tot_ntrks+=v_acxd[i].NbTrk();
if (tot_ntrks<1)
throw ParmError("CxBaselineFitter::CxBaselineFitter(v_data, v_tks) 0 tracks ! tot_ntrks<1 ");
}
int CxBaselineFitter::dofit(string outfilename, string outcheckfilename)
{
size_t NB_ANTENNES=v_acxd[0].getNbAutoCor(); // nombre d'antennes
size_t NB_CXCORS=v_acxd[0].getNbCrossCor();
cout << "======================================================================================"<<endl;
cout << "------- CxBaselineFitter::dofit() Performing baseline/phase fit on the 6 cross-cors "<<" TotNbTracks="<<tot_ntrks<<endl;
POutPersist * pox = NULL ;
if (outcheckfilename.length()>0) {
cout << "... expected cross-cor for fitted params will be saved to file "<<outcheckfilename<<endl;
pox = new POutPersist(outcheckfilename);
}
ofstream ofr(outfilename.c_str());
ofr << "#### Fitted phases and baseline-shifts (CxBaselineFitter::dofit() ) "<<endl
<< "## NumAntenna Phase BaselineShiftX BaselineShiftY BaselineShiftZ (Phase in degree, BaselineShift in meter) "<<endl;
int tot_npoints_fit = 0;
for(size_t i=0; i<v_acxd.size(); i++)
for(size_t j=0; j<v_acxd[i].NbTrk(); j++)
tot_npoints_fit += 2*(v_acxd[i].v_time_data[j].size())*NB_CXCORS;
cout << " Total number of data points for fit="<< tot_npoints_fit<<endl;
GeneralFitData gdata(1, tot_npoints_fit);
for(size_t i=0; i<v_acxd.size(); i++)
for(size_t kcx=0; kcx<NB_CXCORS; kcx++) {
for(size_t j=0; j<v_acxd[i].NbTrk(); j++) {
vector< vector< complex<double> > > & v_cxdata = v_acxd[i].vv_cxdata[j];
vector< vector<double> > & v_cxerr = v_acxd[i].vv_cxerr[j];
for(size_t l=0; l<v_acxd[i].v_time_data[j].size(); l++) {
gdata.AddData1(v_acxd[i].v_time_data[j][l],v_cxdata[kcx][l].real(),v_cxerr[kcx][l]); // Fill x, y and error on y
gdata.AddData1(v_acxd[i].v_time_data[j][l],v_cxdata[kcx][l].imag(),v_cxerr[kcx][l]); // Fill x, y and error on y
}
}
}
My6CxGenXi2B gxi2(v_acxd, v_trks);
GeneralFit mFit(&gxi2);
mFit.SetData(&gdata); // connect data to the fitter , here the data is unused - gxi2 includes its data
mFit.SetMaxStep(1000);
if (NB_ANTENNES != 4)
throw PError("CxBaselineFitter::dofit() NB_ANTENNES != 4 Current version works only for 4 antenna");
// SetParam(int n,double value, double step,double min=1., double max=-1.);
for(size_t i=0; i<(NB_ANTENNES-1); i++) {
char pname[32];
sprintf(pname,"Phase_%d",(int)(i+2));
mFit.SetParam(i,pname,v_acxd[0].v_phase[i],M_PI/360.,0.,2.2*M_PI);
sprintf(pname,"BaselineShift_X_%d",(int)(i+2));
mFit.SetParam(3+3*i,pname,0.,0.01,-0.25,0.25);
sprintf(pname,"BaselineShift_Y_%d",(int)(i+2));
mFit.SetParam(4+3*i,pname,0.,0.01,-0.25,0.25);
sprintf(pname,"BaselineShift_Z_%d",(int)(i+2));
mFit.SetParam(5+3*i,pname,0.,0.01,-0.25,0.25);
}
cout << " Performing the fit ... " << endl;
rcfit = mFit.Fit(); xi2red=-99999.;
cout<< "------ Fit result Reduce_Chisquare = " << mFit.GetChi2Red()<< " nstep="<<mFit.GetNStep() << " rc="<<rcfit<<endl;
mFit.PrintFit();
if (pox) delete pox;
fit_done=true;
return 0;
}
......@@ -85,7 +85,8 @@ public:
double theta_0, phi_0; // Theta, phi angles corresponding
vector< CxBeam > v_acbeams; // the four aut-correlations beams after fit
vector< CxBeam > v_cxbeams; // the six cross correlations after
vector< CxBeam > v_cxbeams; // the six cross correlations after fit
vector<double> v_phase; // fitted phases 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
};
......@@ -160,5 +161,28 @@ public:
vector< CxBeam > v_cxbeams; // the six cross correlations after fit
};
//-----------------------------------------------------------------------
//-------------------------- CxBaselineFitter -------------------------------
//-- Class for fitting baselines and phases using the 6 cross correlations
//-- on a set of tracks
class CxBaselineFitter {
public:
CxBaselineFitter(vector<AcxDataSet> & v_data, vector<TrackSet> & v_tks);
int dofit(string outfilename, string outcheckfilename);
vector<AcxDataSet> & v_acxd;
vector<TrackSet> & v_trks;
size_t tot_ntrks;
bool fit_done;
int rcfit;
double xi2red;
vector<double> v_phases;
vector< Vector3d > v_baselineshits;
};
//-----------------------------------------------------------------------
#endif
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