Skip to content

GitLab

Commit 92a92991 authored by Reza ANSARI's avatar Reza ANSARI
Browse files

Implementation de la possibilite de specifier une shift en coordonnee z des antennes, Reza

parent 335165cf
Hide whitespace changes
Inline Side-by-side
Showing with 99 additions and 21 deletions
trkacxfit.cc
View file @ 92a92991
......@@ -77,6 +77,7 @@ static bool fggaussbeam=true;
// --- if true, perform Beam normalisation
// static bool fgdobeamnormalise=false;
static string outdir_path; // output directory path
static string outfilename; // output filename for auto-correlation fit parameters
static string checkfilename; // output filename for PAON4 and computed auto-cor using fitted parameters
......@@ -86,7 +87,7 @@ static bool fguseAac4Cx=true; //if true, use fitted Amplitudes on autocor
static int prtlevel=0; // print level
static string default_input_path;
static string default_input_path; // default directory for input files
static vector<string> trksetfiles; // datacard files defining the track/data sets
static bool do_cxfit = true ; // if false, perform AutoCor fit only
......@@ -99,7 +100,8 @@ static bool do_baselinefit = false; // if true , perform baseline fit
static bool fg_fixbaseline = false; // if true , do phases fit with fixed baselines
static bool fg_fix_xy = false; // if true , do phases fit with z-baselines free only (fixed x,y-baselines)
static bool do_baselinesimplex = false; // if true , perform baseline determination using Simplex minimisation
//-- shift en Z des antennes
static vector<double> zshift_antennes;
//--- End of list of global parameters
......@@ -119,8 +121,9 @@ int main (int narg, char* arg[])
SophyaInit();
int rcda=decode_args(narg, arg);
if (rcda) return rcda;
TrkFit_SetPrintLevel(prtlevel);
TrkFit_SetZcoordShift(zshift_antennes);
TrkFit_FitLibInfo();
vector<AcxDataSet> v_acxd;
vector<TrackSet> v_trk;
......@@ -133,14 +136,14 @@ int main (int narg, char* arg[])
ACxSetFitter acxfit(acxd, tks);
char fnbuff[32];
sprintf(fnbuff,"ac_%d_",(int)(i+1));
string acofile = fnbuff+outfilename;
string acofile = outdir_path+fnbuff+outfilename;
acxfit.doACfit(acofile);
string acckfile = fnbuff+checkfilename;
string acckfile = outdir_path+fnbuff+checkfilename;
acxfit.saveExpectedAC(acckfile);
if (do_cxfit) { // Fit des cross-cor aussi
sprintf(fnbuff,"cx_%d_",(int)(i+1));
string cxofile = fnbuff+outfilename;
string cxckfile = fnbuff+checkfilename;
string cxofile = outdir_path+fnbuff+outfilename;
string cxckfile = outdir_path+fnbuff+checkfilename;
acxfit.doCxfit(cxofile, fguseAac4Cx, fg_cxB0, fg_phi0only);
acxfit.saveExpectedCx(cxckfile);
}
......@@ -148,44 +151,56 @@ int main (int narg, char* arg[])
v_trk.push_back(tks);
}
CxBaselineFitter cxbfit(v_acxd, v_trk);
if (do_baselinefit || do_baselinesimplex) { // fit simultane des 6 cross-cor
string cxbofile = "cxb6_"+outfilename;
string cxbckfile = "cxb6_"+checkfilename;
CxBaselineFitter cxbfit(v_acxd, v_trk);
int fgfixb=0;
if (fg_fixbaseline) fgfixb=2;
else if (fg_fix_xy) fgfixb=1;
//cxbfit.doSimplexMinimize();
cxbfit.initFitParams();
cxbfit.saveExpectedCx(cxbckfile);
string cxbofile = outdir_path+"cxb6ck_"+outfilename;
string cxbckfile = outdir_path+"cxb6ck_"+checkfilename;
// cxbfit.saveExpectedCx(cxbckfile);
if (do_baselinesimplex) {
cxbfit.doSimplexMinimize(fgfixb);
cxbckfile="cxb6s_"+checkfilename;
cxbckfile=outdir_path+"cxb6s_"+checkfilename;
cxbfit.saveExpectedCx(cxbckfile);
// cxbfit.doCheck();
}
if (do_baselinefit) {
cxbfit.dofit(cxbofile,fgfixb,fg_phi0only);
cxbckfile = "cxb6f_"+checkfilename;
cxbckfile = outdir_path+"cxb6_"+checkfilename;
cxbfit.saveExpectedCx(cxbckfile);
}
}
string summfile=outdir_path+outfilename;
ofstream ofs(summfile.c_str());
{ // Resume resultat fits auto-correlations
cout<<" ---------- Summary Dish diameter / pointing fit from Auto-Correlations ---------- " << endl;
ofs<<" ---------- Summary Dish diameter / pointing fit from Auto-Correlations ---------- " << endl;
for(size_t i=0; i<v_acxd.size(); i++) {
cout << "================================= TrackSet["<<i+1<<"] ================================="<<endl;
ofs << "================================= TrackSet["<<i+1<<"] ================================="<<endl;
v_acxd[i].PrintACFitSummary(cout);
v_acxd[i].PrintACFitSummary(ofs);
}
cout<<"=================================================================================="<<endl;
ofs<<"=================================================================================="<<endl;
}
if (do_cxfit) { // Resume resultat fit cross-cor
cout<<" ---------- Summary phase fitting for cross-correlations ---------- " << endl;
ofs<<" ---------- Summary phase fitting for cross-correlations ---------- " << endl;
for(size_t i=0; i<v_acxd.size(); i++) {
cout << "================================= TrackSet["<<i+1<<"] ================================="<<endl;
ofs << "================================= TrackSet["<<i+1<<"] ================================="<<endl;
v_acxd[i].PrintCxPhaseFitSummary(cout);
v_acxd[i].PrintCxPhaseFitSummary(ofs);
}
cout<<"=================================================================================="<<endl;
ofs<<"=================================================================================="<<endl;
}
cout<<" ---------- Summary 6 cross-correlations baseline fitting ---------- " << endl;
cxbfit.PrintFitSummary(ofs);
Timer tm("trkacxfit");
} // End of try bloc
......@@ -221,15 +236,16 @@ int decode_args(int narg, char* arg[])
if ((narg<2)||fglonghelp) {
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"
<< " options: -inp def_input_path -odp out_path -out OutFileName -ckf CheckFileName \n"
<< " -nocxf -cxfB -phifreq -docx6f -fixb -fixbxy -docx6s \n"
<<" -D dish_diameter -ngb -prt PrintLevel \n"<<endl;
<<" -D dish_diameter -ngb -zshift z1,z2,z3,z4 -prt PrintLevel \n"<<endl;
if (!fglonghelp) {
cout << " trkacxfit -h to get option description " << endl;
return 2;
}
cout << " -inp def_input_path : default directory path for input files (track_Set_J files) \n"
<< " This path can be overwritten in each track_Set_J datacard file \n"
<< " -odp out_path : directory path for output files \n"
<< " -out OutFileName : Output file (text) to save fitted parameters vectors (def trkfit.txt) \n"
<< " -ckf CheckFileName : Output file (PPF) to save expected signals (def chktrkfit.ppf) \n"
<< " Output file names are constructed from OutFileName and CheckFileName prepending \n"
......@@ -244,6 +260,7 @@ int decode_args(int narg, char* arg[])
<< " -docx6s : try simplex minimisation over 6 cross-corr (baseline and phases determination) - 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"
<< " -zshift z1,z2,z3,z4 : define antenna position shift (in cm) along Oz (default:none) \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"
......@@ -252,6 +269,8 @@ int decode_args(int narg, char* arg[])
return 1;
}
default_input_path="";
outdir_path="";
outfilename = "trkfit.txt";
checkfilename = "chktrkfit.ppf";
prtlevel=0;
......@@ -270,6 +289,13 @@ int decode_args(int narg, char* arg[])
if (narg<2) { cout << "trkacxfit/decode_args missing/bad argument, -h for help " << endl; return -1; }
D_dish=atof(arg[2]); arg+=2; narg-=2; lastargs.clear();
}
else if (fbo=="-zshift") { // Antenna shift in z-coord
if (narg<2) { cout << "trkacxfit/decode_args missing/bad argument, -h for help " << endl; return -1; }
double za[4];
sscanf(arg[2],"%lg,%lg,%lg,%lg",za,za+1,za+2,za+3);
zshift_antennes.resize(4); for(size_t ja=0; ja<4; ja++) zshift_antennes[ja]=za[ja];
arg+=2; narg-=2; lastargs.clear();
}
else if (fbo=="-ngb") { // Use Non gaussian beam
fggaussbeam=false; arg++; narg--; lastargs.clear();
}
......@@ -294,10 +320,14 @@ int decode_args(int narg, char* arg[])
else if (fbo=="-docx6s") { // Perform simultaneous fit over 6 cross-corr
do_baselinesimplex=true; arg++; narg--; lastargs.clear();
}
else if (fbo=="-inp") { // output file name
else if (fbo=="-inp") { // default Input path
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();
}
else if (fbo=="-odp") { // output directory path
if (narg<2) { cout << "trkacxfit/decode_args missing/bad argument, -h for help " << endl; return -1; }
outdir_path=arg[2]; arg+=2; narg-=2; lastargs.clear();
}
else if (fbo=="-out") { // output file name
if (narg<2) { cout << "trkacxfit/decode_args missing/bad argument, -h for help " << endl; return -1; }
outfilename=arg[2]; arg+=2; narg-=2; lastargs.clear();
......@@ -322,6 +352,12 @@ 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;
if (zshift_antennes.size()>0) {
cout << " Antenna-Z-coord shift= ";
for(size_t i=0; i<zshift_antennes.size(); i++) cout<<zshift_antennes[i]<<" ; ";
cout<<endl;
}
cout<<" Default input path="<<default_input_path<<" OutputDirectory="<<outdir_path<<endl;
cout << " OutFileName= "<<outfilename<<" CheckFileName= "<<checkfilename<<endl;
cout << " Perform fits on cross-correlations ? " << (do_cxfit?"Yes":"No")
<< " Phase model : Phase(freq)= " << (fg_phi0only?"Phi0":"Phi0+aPhi*(freq-1250)/250")
......
trkfit.cc
View file @ 92a92991
......@@ -27,6 +27,10 @@ void TrkFit_SetPrintLevel(int lev)
return;
}
//--- shift de position en Z pour les 4 antennes
static vector<double> z_coord_shift;
void TrkFit_FitLibInfo()
{
cout << "============================================================================"<<endl;
......@@ -35,10 +39,24 @@ void TrkFit_FitLibInfo()
#else
cout << "============= Classe TkF_Fitter : Fitting with Minuit MnMigrad ============="<<endl;
#endif
if (z_coord_shift.size()>0) {
cout << " Antenna-Z-coord shift= ";
for(size_t i=0; i<z_coord_shift.size(); i++) cout<<z_coord_shift[i]<<" ; ";
cout<<endl;
}
cout << "============================================================================"<<endl;
return;
}
//--- On definit la coordonnees z pour les antennes
void TrkFit_SetZcoordShift(vector<double> & leszs)
{
if (leszs.size()==0) return;
if (leszs.size() != 4) throw ParmError("TrkFit_SetZcoordShift()/ERROR leszs.size() != 4");
z_coord_shift=leszs;
}
//------------------- TrkInputDataSet -------------------------------------
......@@ -703,6 +721,9 @@ int ACxSetFitter::doCxfit(string outfilenamecx, bool useAac, bool fg_B0, bool fg
v_err_Acx[ii][j]=1.; v_err_Bcx[ii][j]=complex<double>(0.,0.);
}
Vector3d baseline=P4Coords::getBaseline(Anum1[ii]+1,Anum2[ii]+1);
if (z_coord_shift.size() > 0) { // Si on a definit un shift des coordonnees z des antennes
baseline += Vector3d(0.,0.,z_coord_shift[Anum2[ii]]-z_coord_shift[Anum1[ii]]);
}
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);
......@@ -957,13 +978,12 @@ CxBaselineFitter::~CxBaselineFitter()
void CxBaselineFitter::initFitParams()
{
// Positions en z des 3 antennes 2,3,4 / 1 , en metres
double z_antennes[3]={-0.10,+0.05,-0.02};
//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)
throw PError("CxBaselineFitter::initFitParams() NB_ANTENNES != 4 Current version works only for 4 antenna");
v_phase.resize(v_acxd[0].getNbAutoCor()-1);
v_phi_0.resize(v_acxd[0].getNbAutoCor()-1);
v_err_phi_0.resize(NB_ANTENNES-1);
v_a_phi.resize(v_acxd[0].getNbAutoCor()-1);
......@@ -974,6 +994,7 @@ void CxBaselineFitter::initFitParams()
v_phi_0[i]=v_acxd[0].v_phi_0[i]; v_err_phi_0[i]=0.;
v_a_phi[i]=v_acxd[0].v_a_phi[i]; v_err_a_phi[i]=0.;
v_baselineshits[i]=Vector3d(0.,0.,0.);
if (z_coord_shift.size()>0) v_baselineshits[i]=Vector3d(0.,0.,z_coord_shift[i+1]-z_coord_shift[0]);
v_err_baselineshits[i]=Vector3d(0.,0.,0.);
bestfitparam[2*i]=v_phi_0[i];
err_bestfitparam[2*i]=0.;
......@@ -982,6 +1003,7 @@ void CxBaselineFitter::initFitParams()
for(size_t j=0; j<3; j++) {
bestfitparam[3*(i+1)+j]=err_bestfitparam[3*(i+1)+j]=0.;
}
if (z_coord_shift.size()>0) bestfitparam[3*(i+1)+2]=z_coord_shift[i+1]-z_coord_shift[0];
}
//DBG cout << " *DBG* DONE **** CxBaselineFitter::initFitParams()"<<endl;
......@@ -1050,7 +1072,7 @@ int CxBaselineFitter::dofit(string outfilename, int fgfixbaseline, bool fgphi0on
}
}
int ndataused=0;
double xi2redini=gxi2.getXi2(bestfitparam, ndataused);
xi2redini=gxi2.getXi2(bestfitparam, ndataused);
cout << " Xi2Red value for the initial parameters guess= " <<xi2redini/(double)ndataused<<endl;
cout << " Performing the fit (tot_npoints_fit= "<<tot_npoints_fit<<" ?= (npoints2="<<npoints2<<") ..."<< endl;
rcfit = mFit.doFit(); xi2red=-99999.;
......@@ -1080,6 +1102,7 @@ int CxBaselineFitter::dofit(string outfilename, int fgfixbaseline, bool fgphi0on
double phase=gxi2.getPhase4Freq(v_phi_0[i],v_a_phi[i],1300.);
while (phase<0.) phase += 2.*M_PI;
while (phase>2.*M_PI) phase -= 2.*M_PI;
v_phase[i]=phase;
cout <<"---Ant["<<i+2<<" Phase(@1300MHz)= "<<setw(10)<<Angle(phase).ToDegree()<<" phi_0= "<<setw(10)
<<Angle(v_phi_0[i]).ToDegree()<<" +/- "<<setw(10)<<Angle(v_err_phi_0[i]).ToDegree()<<" deg."
<<" a_phi= "<<setw(8)<<Angle(v_a_phi[i]).ToDegree()<<" +/- "<<setw(10)
......@@ -1090,6 +1113,21 @@ int CxBaselineFitter::dofit(string outfilename, int fgfixbaseline, bool fgphi0on
return 0;
}
ostream & CxBaselineFitter::PrintFitSummary(ostream & os)
{
os<<"--------------------------------------------------------------------------------------------"<<endl;
os<< "---CxBaselineFitter result Reduce_Chisquare = "<<xi2red<<" rc="<<rcfit<<" Xi2RedInitial="<<xi2redini<<endl;
for(size_t i=0; i<v_acxd[0].getNbAutoCor()-1; i++) {
os <<"---Ant["<<i+2<<" Phase(@1300MHz)= "<<setw(10)<<Angle(v_phase[i]).ToDegree()<<" phi_0= "<<setw(10)
<<Angle(v_phi_0[i]).ToDegree()<<" +/- "<<setw(10)<<Angle(v_err_phi_0[i]).ToDegree()<<" deg."
<<" a_phi= "<<setw(8)<<Angle(v_a_phi[i]).ToDegree()<<" +/- "<<setw(10)
<<Angle(v_err_a_phi[i]).ToDegree()<<" deg/250 MHz"<<endl;
os <<" Baselineshift["<<i+2<<" = "<<v_baselineshits[i] << " +/- " << v_err_baselineshits[i] << " m."<<endl;
}
os<<"--------------------------------------------------------------------------------------------"<<endl;
return os;
}
int CxBaselineFitter::doSimplexMinimize(int fg_fixebaseline)
{
size_t NB_ANTENNES=v_acxd[0].getNbAutoCor(); // nombre d'antennes
......@@ -1122,7 +1160,7 @@ int CxBaselineFitter::doSimplexMinimize(int fg_fixebaseline)
}
cout << " Initial Point: "<<P0.Transpose()<<endl;
cout << " Initial Step: "<<step.Transpose()<<endl;
cout << " Initial Xi2= " << mzfunc.Value(P0.Data())<<endl;
cout << " ==> Initial Xi2= " << mzfunc.Value(P0.Data())<<endl;
simplex.SetInitialPoint(P0);
simplex.SetInitialStep(step);
......
trkfit.h
View file @ 92a92991
......@@ -22,6 +22,7 @@ using namespace SOPHYA;
void TrkFit_SetPrintLevel(int lev=0);
void TrkFit_FitLibInfo();
void TrkFit_SetZcoordShift(vector<double> & leszs);
//-------------------------------------------------------------------------
//------------------- TrkInputDataSet -------------------------------------
......@@ -214,6 +215,8 @@ public:
int saveExpectedCx(string outcheckfilename); // outfilename PPF (or fits ?) file with the expected signals
ostream & PrintFitSummary(ostream & os);
vector<AcxDataSet> & v_acxd;
vector<TrackSet> & v_trks;
size_t tot_ntrks;
......@@ -221,7 +224,8 @@ public:
bool fit_done;
bool simplex_done;
int rcfit;
double xi2red;
double xi2redini, xi2red;
vector<double> v_phase;
vector<double> v_phi_0;
vector<double> v_err_phi_0;
vector<double> v_a_phi;
......
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