Skip to content

GitLab

Commit d0c29e2c authored by Reza ANSARI's avatar Reza ANSARI
Browse files

Merge branch 'master' of gitlab.in2p3.fr:baoradio/AnaPAON4 

  Merge avec les modifs d'Olivier sur le programme de fit - j'ai juste modifie
  des scripts   Reza
parents 4088e8ee b99382b3
Hide whitespace changes
Inline Side-by-side
Showing with 207 additions and 27 deletions
Python/viewtfmap.py
View file @ d0c29e2c
......@@ -70,6 +70,7 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
if timin ==0. :
timin=times[0]/3600.
timax = max(times)/3600.
print ' timin/max = '+str(timin)+', ' +str(timax)
if ramod :
timin = min(ras)
timax= max(ras)
......@@ -153,7 +154,8 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
print vmax,np.asarray(img[300:900,]).max()
if (not noplo) :
im = ax1.imshow(img,aspect='auto' ,extent=(timin,timax,frmax,frmin), interpolation='none',vmin=vmin,vmax=vmax)
im = ax1.imshow(img,aspect='auto' , interpolation='none',vmin=vmin,vmax=vmax)
#m = ax1.imshow(img,aspect='auto' ,extent=(timin,timax,frmax,frmin), interpolation='none',vmin=vmin,vmax=vmax)
if (yrcm) :
im.set_cmap('YlOrBr_r')
cbar = plt.colorbar( im ) #,fraction=0.021,pad=.05)
......@@ -180,13 +182,13 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
ax1.annotate(nmras[plras.index(rap)],(rap,frmin),horizontalalignment='left',verticalalignment='bottom',rotation=30)
ax1.yaxis.grid(which="major", color='gray', linestyle='-', linewidth=.5)
ax1.xaxis.grid(which="major", color='gray', linestyle='-', linewidth=.5)
ax2=ax1.twinx()
ax2.yaxis.set_view_interval(len(freqs),0.,ignore=True)
ax2.yaxis.set_view_interval(frmax,frmin ,ignore=True)
ax2.set_ylabel("Freq. bin",rotation=270.,y=0.9)
if not ramod :
ax3=ax1.twiny()
ax3.xaxis.set_view_interval(0.,len(times))
ax3.xaxis.set_view_interval( timin,timax ,ignore=True)
ax3.set_xlabel("time bin",x=0.85)
if save:
......@@ -254,7 +256,7 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
bnext.on_clicked(cm_next)
bprev.on_clicked(cm_prev)
plt.tight_layout()
#plt.tight_layout()
plt.subplots_adjust(bottom=0.15)
plt.show()
......
Satellites/trk2dt.cc
View file @ d0c29e2c
......@@ -89,6 +89,11 @@ int main(int narg, const char* arg[])
DataTableRow row = dt.EmptyRow();
ifstream is(infile.c_str());
if (is.fail()){
cout<<" Problem opening file !! "<<endl;
throw IOExc(" trk2dt : bad file : "+infile+" - theck it !!");
}
string sline;
while (!is.eof()) {
sline = "";
......@@ -141,7 +146,7 @@ int main(int narg, const char* arg[])
}
}
catch (PException& exc) {
cerr << " rdvisip4.cc catched PException " << exc.Msg() << endl;
cerr << " trk2dt.cc catched PException " << exc.Msg() << endl;
rc = 77;
}
catch (std::exception& sex) {
......
trkfit.cc
View file @ d0c29e2c
......@@ -7,6 +7,8 @@
#include <iomanip>
#include "pexceptions.h"
#include "trkfit.h"
#include "datacards.h"
#include "array.h"
......@@ -17,6 +19,9 @@
#include "gcxfitbaseline.h"
#include "p4autils.h"
#include "fitsioserver.h"
using namespace std;
using namespace SOPHYA;
//------------------- Print Level for this file --------------------------
......@@ -90,6 +95,10 @@ static int decode_trkcard(string const& key, string const& toks)
char flnmdata[256], flnmtrk[256];
char sflags[64];
double ts,te,freq;
flnmdata[0]='\0';
flnmtrk[0] ='\0';
sflags[0] ='\0';
ts=te=freq=-9999.;
sscanf(toks.c_str(),"%s %lg,%lg %lg %s %s",flnmdata,&ts,&te,&freq,flnmtrk,sflags);
dataflnm_p_->push_back(flnmdata);
......@@ -307,7 +316,7 @@ size_t AcxDataSet::ReadData(TrkInputDataSet & tkds)
vector< complex<double> > & vcx = v_vcx[ii];
v_cxdata[ii].push_back(vcx[k]);
double acx=std::abs(vcx[k]);
v_cxerr[ii].push_back(0.1*sqrt(acx));
v_cxerr[ii].push_back(0.1*sqrt(acx)); // calcul d'erreur, a affiner
if (acx<v_min_cxdata[j][ii]) v_min_cxdata[j][ii]=acx;
if (acx>v_max_cxdata[j][ii]) v_max_cxdata[j][ii]=acx;
}
......@@ -649,10 +658,21 @@ int ACxSetFitter::doACfit(string outfilename)
int ACxSetFitter::saveExpectedAC(string outcheckfilename)
{
if (outcheckfilename.length()<1) return 1;
cout << "-----ACxSetFitter::saveExpectedAC() : computing expected signal for fitted params , will be saved to file "
<<outcheckfilename<<endl;
POutPersist pos(outcheckfilename);
vector <string> data_names;
FitsInOutFile *fos=NULL;
POutPersist *pos=NULL;
size_t ll = outcheckfilename.length();
size_t pp = outcheckfilename.rfind('.');
if ((pp<ll)&&(outcheckfilename.substr(pp)==".fits")) {
cout << "----- ACxSetFitter::saveExpectedAC : Saving to FITS file "<<outcheckfilename<<endl;
fos=new FitsInOutFile (outcheckfilename,FitsInOutFile::Fits_Create);
} else {
pos = new POutPersist (outcheckfilename);
}
size_t NB_ANTENNES = acxd_.getNbAutoCor();
size_t NTRK = acxd_.NbTrk();
......@@ -671,19 +691,56 @@ int ACxSetFitter::saveExpectedAC(string outcheckfilename)
for(size_t j=0; j<NTRK; j++) {
double A = v_A[ii][j];
double B = v_B[ii][j];
Vector signal = macs.getDataSignal(j);
sprintf(oname,"ac_%d_%d",(int)ii+1,(int)j+1);
pos << PPFNameTag(oname)<<signal;
Vector expsignal = macs.getExpectedSignal(j, Ddishfit, thetafit, phifit, A, B);
sprintf(oname,"simac_%d_%d",(int)ii+1,(int)j+1);
pos << PPFNameTag(oname)<<expsignal;
if (ii==0) {
Vector tmvec = macs.getTimeVec(j);
sprintf(oname,"tim_%d",(int)j+1);
pos << PPFNameTag(oname)<<tmvec;
if (pos!=NULL) {
*pos << PPFNameTag(oname)<<tmvec;
}
if (fos!=NULL) {
data_names.push_back(oname);
tmvec.Info()["NOM_OBJET"]=oname;
*fos<<tmvec;
}
}
Vector signal = macs.getDataSignal(j);
sprintf(oname,"ac_%d_%d",(int)ii+1,(int)j+1);
if (pos!=NULL) {
*pos << PPFNameTag(oname)<<signal;
}
if (fos!=NULL) {
data_names.push_back(oname);
signal.Info()["NOM_OBJET"]=oname;
*fos<<signal;
}
Vector expsignal = macs.getExpectedSignal(j, Ddishfit, thetafit, phifit, A, B);
sprintf(oname,"simac_%d_%d",(int)ii+1,(int)j+1);
if (pos!=NULL) {
*pos << PPFNameTag(oname)<<expsignal;
}
if (fos!=NULL) {
expsignal.Info()["NOM_OBJET"]=oname;
data_names.push_back(oname);
*fos<<expsignal;
}
}
}
if (pos!=NULL) delete(pos);
if (fos!=NULL) {
//*fos<<data_names;
//TVector <string> tata(data_names);
//*fos<<tata ;
cout<<" number of extenstions "<<data_names.size()<<endl;
for (int k=0 ; k<data_names.size() ; k++){
cout<<" no "<<k<<" ->"<<data_names[k]<<endl;
}
delete(fos);
}
return 0;
}
......@@ -919,7 +976,19 @@ int ACxSetFitter::doCxfit(string outfilenamecx, bool useAac, bool fg_B0, bool fg
int ACxSetFitter::saveExpectedCx(string outcheckfilename)
{
cout << "ACxSetFitter::saveExpectedCx() saving expected cross-cor (and visi-data) to file "<<outcheckfilename<<endl;
POutPersist pox(outcheckfilename);
size_t ll = outcheckfilename.length();
size_t pp = outcheckfilename.rfind('.');
FitsInOutFile *fox=NULL;
POutPersist *pox=NULL;
vector <string> data_names;
if ((pp<ll)&&(outcheckfilename.substr(pp)==".fits")) {
cout << "----- ACxSetFitter::saveExpectedCx : Saving to FITS file "<<outcheckfilename<<endl;
fox=new FitsInOutFile (outcheckfilename,FitsInOutFile::Fits_Create);
} else {
pox = new POutPersist (outcheckfilename);
}
size_t NB_CXCORS=acxd_.getNbCrossCor();
size_t NTRK = acxd_.NbTrk();
......@@ -930,21 +999,69 @@ int ACxSetFitter::saveExpectedCx(string outcheckfilename)
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);
if (ii==0) {
Vector tmvec = cxsig.getTimeVec(j);
sprintf(oname,"tim_%d",(int)j+1);
if (pox!=NULL)
*pox << PPFNameTag(oname)<<tmvec;
if (fox!=NULL){
data_names.push_back(oname);
tmvec.Info()["NOM_OBJET"]=oname;
*fox << tmvec;
}
}
TVector < complex<double> > signal = cxsig.getDataSignal(j);
sprintf(oname,"cx_%d_%d",(int)ii+1,(int)j+1);
pox << PPFNameTag(oname)<<signal;
if (pox!=NULL)
*pox << PPFNameTag(oname)<<signal;
if (fox!=NULL){
sprintf(oname,"cx_%d_%d_R",(int)ii+1,(int)j+1);
data_names.push_back(oname);
Vector tmpv=real(signal);
tmpv.Info()["NOM_OBJET"]=oname;
*fox << tmpv ;
sprintf(oname,"cx_%d_%d_I",(int)ii+1,(int)j+1);
data_names.push_back(oname);
tmpv=imag(signal);
tmpv.Info()["NOM_OBJET"]=oname;
*fox << tmpv ;
}
//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) {
Vector tmvec = cxsig.getTimeVec(j);
sprintf(oname,"tim_%d",(int)j+1);
pox << PPFNameTag(oname)<<tmvec;
if (pox!=NULL)
*pox << PPFNameTag(oname)<<expsignal;
if (fox!=NULL){
sprintf(oname,"simcx_%d_%d_R",(int)ii+1,(int)j+1);
data_names.push_back(oname);
Vector tmpv = real(expsignal);
tmpv.Info()["NOM_OBJET"]=oname;
*fox << tmpv ;
sprintf(oname,"simcx_%d_%d_I",(int)ii+1,(int)j+1);
data_names.push_back(oname);
tmpv = imag(expsignal);
tmpv.Info()["NOM_OBJET"]=oname;
*fox << tmpv ;
}
}
}
if (pox!=NULL) delete(pox);
if (fox!=NULL) {
cout<<" number of extensions "<<data_names.size()<<endl;
for (int k=0 ; k<data_names.size() ; k++){
cout<<" no "<<k<<" ->"<<data_names[k]<<endl;
}
delete(fox);
}
return 0;
}
......@@ -1272,7 +1389,18 @@ int CxBaselineFitter::doCheck(int fg_fixebaseline)
int CxBaselineFitter::saveExpectedCx(string outcheckfilename)
{
cout << "CxBaselineFitter::saveExpectedCx() saving expected cross-cor (and visi-data) to file "<<outcheckfilename<<endl;
POutPersist pox(outcheckfilename);
FitsInOutFile *fox=NULL;
POutPersist *pox=NULL;
vector <string> data_names;
size_t ll = outcheckfilename.length();
size_t pp = outcheckfilename.rfind('.');
if ((pp<ll)&&(outcheckfilename.substr(pp)==".fits")) {
cout << "----- CxBaselineFitter::saveExpectedCx : Saving to FITS file "<<outcheckfilename<<endl;
fox = new FitsInOutFile (outcheckfilename,FitsInOutFile::Fits_Create);
} else {
pox = new POutPersist (outcheckfilename);
}
size_t NB_CXCORS=v_acxd[0].getNbCrossCor();
char oname[48];
......@@ -1283,16 +1411,61 @@ int CxBaselineFitter::saveExpectedCx(string outcheckfilename)
for(size_t j=0; j<v_acxd[i].NbTrk(); j++) {
Vector tmvec = cxsigb.getTimeVec(i,j);
sprintf(oname,"tim_%d_%d",(int)j+1,(int)i+1);
pox << PPFNameTag(oname)<<tmvec;
if (pox!=NULL)
*pox << PPFNameTag(oname)<<tmvec;
if (fox!=NULL){
data_names.push_back(oname);
tmvec.Info()["NOM_OBJET"]=oname;
*fox << tmvec;
}
for(size_t kcx=0; kcx<NB_CXCORS; kcx++) {
TVector< complex<double> > signal = cxsigb.getDataSignal(i,j,kcx);
sprintf(oname,"cx_%d_%d_%d",(int)kcx+1,(int)j+1,(int)i+1);
pox << PPFNameTag(oname)<<signal;
if (pox!=NULL){
*pox << PPFNameTag(oname)<<signal;
}
if (fox!=NULL){
sprintf(oname,"cx_%d_%d_%d_R",(int)kcx+1,(int)j+1,(int)i+1);
data_names.push_back(oname);
Vector vtmp=real(signal);
vtmp.Info()["NOM_OBJET"]=oname;
*fox << vtmp ;
sprintf(oname,"cx_%d_%d_%d_I",(int)kcx+1,(int)j+1,(int)i+1);
data_names.push_back(oname);
vtmp=imag(signal);
vtmp.Info()["NOM_OBJET"]=oname;
*fox << imag(signal);
}
TVector< complex<double> > expsignal = cxsigb.getExpectedSignal(i,j,kcx,bestfitparam);
sprintf(oname,"simcx_%d_%d_%d",(int)kcx+1,(int)j+1,(int)i+1);
pox << PPFNameTag(oname)<<expsignal;
if (pox!=NULL){
*pox << PPFNameTag(oname)<<expsignal;
}
if (fox!=NULL){
sprintf(oname,"simcx_%d_%d_%d_R",(int)kcx+1,(int)j+1,(int)i+1);
data_names.push_back(oname);
Vector vtmp=real(expsignal);
vtmp.Info()["NOM_OBJET"]=oname;
*fox << vtmp;
sprintf(oname,"simcx_%d_%d_%d_I",(int)kcx+1,(int)j+1,(int)i+1);
data_names.push_back(oname);
vtmp=imag(expsignal);
vtmp.Info()["NOM_OBJET"]=oname;
*fox << vtmp ;
}
}
}
}
if (pox!=NULL) delete(pox);
if (fox!=NULL) {
cout<<" number of extenstions "<<data_names.size()<<endl;
for (int k=0 ; k<data_names.size() ; k++){
cout<<" no "<<k<<" ->"<<data_names[k]<<endl;
}
delete(fox);
}
return 0;
}
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