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

Added computation of frequency dependent gain g(nu) in rdvisip4.cc , and use...

Added computation of frequency dependent gain g(nu) in rdvisip4.cc , and use of g(nu) in visi2ntac.cc, as well as a second frequency band in the N-tuple, Reza 13/04/2015
parent 7557c090
......@@ -9,6 +9,10 @@ List of files:
produced by mfac, merge the different frequency bands and set of visibilities
and reorder frequencies
- rdvisip4.cc : a simple main program to read visibilitiy files and produce a single
mean visibility matrix
mean visibility matrix
Updated April 2015: computes also a frequency dependent gain gain(nu) using a gain
function developed by Christophe (cmv) for cluster data analysis
- visi2ntac.cc : a simple main program to read visibilitiy files and produce a NTuple
with total power as a function of time for the 8 PAON4 auto-correlations
Updated April 2015: reads the gain g(nu) file produced by rdvisip4.cc , corrects the
auto-correlation signals before integrating in the specified bands
......@@ -39,6 +39,13 @@
// include lecteur de fichiers visibilites
#include "visip4reader.h"
// ---- compute gain for each channel from visibilies and save it to file
void computeSaveGain(TMatrix< complex<r_4> > & vismtx_mean, string const & gainfile);
// ----- adaptation de la fonction de calcul de gain, ecrit par cmv pour l'analyse Amas (Aout 2013)
Matrix computeGain2(Matrix & acs, int SZW=96,float frackeep=0.75);
// ----- fonction pour normaliser les gains
Vector normalizeGain2(Matrix& gain);
//--------------------------- Fonctions de ce fichier -------------------
int Usage(bool fgshort=true);
// int DecodeArgs(int narg, char* arg[]);
......@@ -48,10 +55,11 @@ int Usage(bool fgea)
{
cout << " --- rdvisip4.cc : Read PPF files produced by mfacq time-frequency\n" << endl;
if (fgea) cout << " rdvisip4: Missing or wrong argument ! " << endl;
cout << " Usage: rdvisip4 InPathBAO5 InPathBAO6 Imin,Imax,step OutPPFFile [PrintLev=1] \n"
cout << " Usage: rdvisip4 InPathBAO5 InPathBAO6 Imin,Imax,step OutPPFFile GainPPFFile [PrintLev=1] \n"
<< " InPathBAO5/6: path for input vismtx_J_iii.ppf on bao5/6 \n"
<< " Imin,Imax,step : read files for Imin<=iii<=Imax iii+=step \n"
<< " OutPPFFile: Output PPF file name \n"<<endl;
<< " OutPPFFile: Output PPF file name (with the average visibility matrix) \n"
<< " GainPPFFile: PPF file name containing the computed gain g(nu) \n"<<endl;
/*
if (fgshort) {
cout << " rdvisip4 -h for detailed instructions" << endl;
......@@ -66,14 +74,15 @@ int Usage(bool fgea)
int main(int narg, char* arg[])
{
if ((narg>1)&&(strcmp(arg[1],"-h")==0)) return Usage(false);
if (narg<5) return Usage(true);
if (narg<6) return Usage(true);
string path5 = arg[1];
string path6 = arg[2];
int Imin,Imax,Istep;
sscanf(arg[3],"%d,%d,%d",&Imin,&Imax,&Istep);
string outfile=arg[4];
string gainfile=arg[5];
int prtlev=1;
if (narg>5) prtlev=atoi(arg[5]);
if (narg>6) prtlev=atoi(arg[6]);
cout << " rdvisip4/Info: Path BAO5:"<<path5<<" BAO6:"<<path6<<"\n"
<< " Imin,max,step="<<Imin<<","<<Imax<<","<<Istep<<" OutFile:"<<outfile
......@@ -108,6 +117,8 @@ int main(int narg, char* arg[])
cout << " rdvisip4/Info: count="<<cnt<<" visimtx read "<<endl;
if (cnt>0) {
vismtx_mean /= complex<r_4>((r_4)cnt, (r_4)0.);
cout << " rdvisip4/Info: computing gains ..."<<endl;
computeSaveGain(vismtx_mean, gainfile);
cout<<" rdvisip4/Info: Saving vismtx_mean to PPF file "<<outfile<<endl;
POutPersist po(outfile);
po<<vismtx_mean;
......@@ -135,4 +146,91 @@ int main(int narg, char* arg[])
}
//--------------------- Fonctions pour le calcul de gain -----------------------
//----- include specifique pour computeGain2()
#include "nbtrixx.h"
#include "cspline.h"
/*--Fonction--*/
void computeSaveGain(TMatrix< complex<r_4> > & vismtx_mean, string const & gainfile)
{
sa_size_t KVAC[8]={0,8,15,21,26,30,33,35};
Matrix acs(8, vismtx_mean.NCols());
for(sa_size_t r=0; r<8; r++)
for(sa_size_t c=0; c<acs.NCols(); c++) acs(r,c)=vismtx_mean(KVAC[r],c).real();
int SZW=96; float frackeep=0.75;
cout << " computeSaveGain/Info: computing gains with SZW="<<SZW<<" frackeep="<<frackeep<<" ..."<<endl;
Matrix gains = computeGain2(acs, SZW, frackeep);
Vector gn = normalizeGain2(gains);
cout << " computeSaveGain/Info: saving gains, gn to "<<gainfile<<endl;
POutPersist po(gainfile);
po<<PPFNameTag("gains")<<gains;
po<<PPFNameTag("gn")<<gn;
return;
}
// ----- adaptation de la fonction de calcul de gain, ecrit par cmv pour l'analyse Amas (Aout 2013)
/* --- fonction de calcul de gain avec fit de spline qui flag tout seul
les mauvais points du spectre - calcule la moyenne par fenetre (en
utilisant une mediane) et fit de spline sur les un point / fenetre
base sur la fonction rawgain2 de cmv (en date du gain (fonction de normalisation)
des spectres ON et OFF separes
*/
/*--Fonction--*/
Matrix computeGain2(Matrix& acs, int SZW,float frackeep)
{
if(SZW%2) SZW++;
int nrows = acs.NRows(), ncols = acs.NCols();
double *val = new double[2*SZW+2];
Matrix gain(nrows,ncols); gain = 0.;
for(int r=0;r<nrows;r++) {
TVector<r_8> DVon(ncols); DVon = 0.;
for(int c=1; c<ncols-1; c++) DVon(c) = fabs(acs(r,c) - 0.5*(acs(r,c-1)+acs(r,c+1)));
vector<double> vsplx, vsply;
for(int c=0; c<ncols; c+=SZW/4) {
// on enleve d'abord les RFI tres fortes
int ng = 0;
for(int cc=c;cc<c+SZW&&cc<ncols;cc++) {val[ng] = DVon(cc); ng++;}
if(ng<2) continue;
TabSort(ng,val);
double dvcut = val[int(ng*frackeep)];
//
ng = 0;
double x = 0.;
for(int cc=c;cc<c+SZW&&cc<ncols;cc++) if(DVon(cc)<dvcut) {val[ng] = acs(r,cc); x+=cc; ng++;}
if(ng!=0) x /= (double)ng;
////cout<<"rawgain2: r="<<r<<" c="<<c<<" dvcut="<<dvcut<<" ng="<<ng<<" x="<<x<<" v="<<val[ng/2]<<endl;
if(ng<2) continue;
TabSort(ng,val);
vsplx.push_back(x);
vsply.push_back(val[ng/2]);
}
TVector<r_8> vvx(vsplx), vvy(vsply);
CSpline csp(vvx.Size(),vvx.Data(),vvy.Data());
csp.ComputeCSpline();
for(int c=0;c<ncols;c++) gain(r,c) = csp.CSplineInt(c);
}
delete [] val;
}
/*--Fonction--*/
Vector normalizeGain2(Matrix& gain)
{
char buff[16];
Vector gnorm(gain.NRows());
for(sa_size_t r=0;r<gain.NRows();r++) {
double gn=gain.Row(r).Sum();
gain.Row(r) /= gn;
sprintf(buff,"GNORM%d",(int)r);
gain.Info()[buff]=gn;
gnorm(r)=gn;
size_t badcnt=0;
#define LOWGTHR 1.e-9
for(sa_size_t c=0; c<gain.NCols(); c++)
if (gain(r,c)<LOWGTHR) {
badcnt++; gain(r,c)=LOWGTHR;
}
cout << " normalizeGain2()/Info GNorm["<<r<<"]= "<<gn<<" (badcount="<<badcnt<<")"<<endl;
}
return gnorm;
}
......@@ -48,11 +48,12 @@ int Usage(bool fgea)
{
cout << " --- visi2ntac.cc : Read PPF files produced by mfacq time-frequency\n" << endl;
if (fgea) cout << " visi2ntac: Missing or wrong argument ! " << endl;
cout << " Usage: visi2ntac InPathBAO5 InPathBAO6 Imin,Imax OutPPFFile [Jf1,Jf2] [DeltaIAvg] [PrtLev=0] \n"
cout << " Usage: visi2ntac InPathBAO5 InPathBAO6 Imin,Imax GainPPFFile OutPPFFile Jf1A,Jf2A JF1B,Jf2B DeltaIAvg [PrtLev=0] \n"
<< " InPathBAO5/6: path for input vismtx_J_iii.ppf on bao5/6 \n"
<< " Imin,Imax : read files for Imin<=iii<=Imax iii+=step \n"
<< " GainPPFFile: Input gains PPF file name \n"
<< " OutPPFFile: Output PPF file name \n"
<< " Jf1,Jf2: frequency index (Jf1<=f<=Jf2) range for power integration, def=768,2304 \n"
<< " Jf1A/B,Jf2A/B : frequency index (Jf1<=f<=Jf2) range for power integration, A:1300,1460 B:2130,2290 \n"
<< " DeltaIAvg: compute average power every DeltaI input vismtx files, def=1 \n"<<endl;
/*
if (fgshort) {
......@@ -68,29 +69,44 @@ int Usage(bool fgea)
int main(int narg, char* arg[])
{
if ((narg>1)&&(strcmp(arg[1],"-h")==0)) return Usage(false);
if (narg<5) return Usage(true);
if (narg<9) return Usage(true);
string path5 = arg[1];
string path6 = arg[2];
int Imin,Imax,Istep=1;
sscanf(arg[3],"%d,%d",&Imin,&Imax);
Istep=1;
string outfile=arg[4];
string gainfile=arg[4];
string outfile=arg[5];
// frequency range definition
sa_size_t JFmin=768,JFmax=2304;
sa_size_t JFmin=1300,JFmax=1460;
int jf1,jf2;
if ((narg>5)&&(strcmp(arg[5],".")!=0)) sscanf(arg[5],"%d,%d",&jf1,&jf2);
if ((narg>6)&&(strcmp(arg[6],"-")!=0)) sscanf(arg[6],"%d,%d",&jf1,&jf2);
JFmin=jf1; JFmax=jf2;
sa_size_t JFminB=2130,JFmaxB=2290;
if ((narg>7)&&(strcmp(arg[7],"-")!=0)) sscanf(arg[7],"%d,%d",&jf1,&jf2);
JFminB=jf1; JFmaxB=jf2;
// time averaging interval definition, in term of visibility files
int deltaIavg=1;
if (narg>6) deltaIavg=atoi(arg[6]);
int deltaIavg=10;
if (narg>8) deltaIavg=atoi(arg[8]);
if (deltaIavg<1) deltaIavg=1;
int prtlev=1;
if (narg>7) prtlev=atoi(arg[7]);
if (narg>9) prtlev=atoi(arg[9]);
cout << " visi2ntac/Info: Path BAO5:"<<path5<<" BAO6:"<<path6<<"\n"
<< " Imin,max,step="<<Imin<<","<<Imax<<","<<Istep<<" OutFile:"<<outfile<<"\n"
<< JFmin<<" <=NumFreq<= "<<JFmax<<" DeltaIAvg="<<deltaIavg<<" PrtLev="<<prtlev<<endl;
<< JFmin<<" <=NumFreq<= "<<JFmax<<" "<<JFminB<<" <=NumFreqB<= "<<JFmaxB
<<" DeltaIAvg="<<deltaIavg<<" PrtLev="<<prtlev<<endl;
cout << " visi2ntac/Info: reading input gain file"<<gainfile<<endl;
Matrix gains;
Vector gn;
{
PInPersist pin(gainfile);
pin>>PPFNameTag("gains")>>gains;
pin>>PPFNameTag("gn")>>gn;
for (int p=0; p<8; p++) cout<<"gn["<<p<<"]="<<gn(p)<<" ";
cout<<endl;
}
// Bande de frequence de 2 MHz autour du 21 cm (1420 MHz -> JF=2792)
sa_size_t JFmin21=2776,JFmax21=2808; // +/- 1 MHz autour de 1420 MHz
sa_size_t JFmin21_5=2710,JFmax21_5=2874; // +/- 5 MHz autour de 1420 MHz
......@@ -107,18 +123,20 @@ int main(int narg, char* arg[])
paths.push_back(path5);
paths.push_back(path6);
const char* ntnames[26]={"k","time",
const char* ntnames[34]={"k","time",
"p1","p2","p3","p4","p5","p6","p7","p8",
"p1HI","p2HI","p3HI","p4HI","p5HI","p6HI","p7HI","p8HI",
"p1HI5","p2HI5","p3HI5","p4HI5","p5HI5","p6HI5","p7HI5","p8HI5"};
"p1HI5","p2HI5","p3HI5","p4HI5","p5HI5","p6HI5","p7HI5","p8HI5",
"p1B","p2B","p3B","p4B","p5B","p6B","p7B","p8B" };
NTuple nt(26,ntnames);
r_8 xnt[30];
NTuple nt(34,ntnames);
r_8 xnt[40];
// Numero de ligne des auto-correlations dans la matrice des visibilites
sa_size_t KVAC[8]={0,8,15,21,26,30,33,35};
Range freqs(JFmin,JFmax);
Range freqs21(JFmin21,JFmax21);
Range freqs21_5(JFmin21_5,JFmax21_5);
Range freqsB(JFminB,JFmaxB);
VisiP4Reader vreader(paths, Imin,Imax,Istep,true);
vreader.setPrintLevel(prtlev);
......@@ -138,10 +156,14 @@ int main(int narg, char* arg[])
for(int k=0; k<24; k++) xnt[2+k]=0.;
}
for(int k=0; k<8; k++) {
TVector< complex<r_4> > vac = vismtx.Row(KVAC[k]);
xnt[2+k] += vac(freqs).Sum().real(); // integration en bande large
xnt[2+8+k] += vac(freqs21).Sum().real(); // integration a +/- 1 MHz @ 1420 MHz
xnt[2+16+k] += vac(freqs21_5).Sum().real(); // integration a +/- 5 MHz @ 1420 MHz
TVector< complex<r_4> > vacz = vismtx.Row(KVAC[k]);
Vector vac(vismtx.NCols());
// correct for gain(nu)
for(sa_size_t c=0; c<vac.Size(); c++) vac(c)=vacz(c).real()/gains(k,c);
xnt[2+k] += vac(freqs).Sum(); // integration en bande large
xnt[2+8+k] += vac(freqs21).Sum(); // integration a +/- 1 MHz @ 1420 MHz
xnt[2+16+k] += vac(freqs21_5).Sum(); // integration a +/- 5 MHz @ 1420 MHz
xnt[2+24+k] += vac(freqsB).Sum(); // integration en bande large B
}
I++;
if (I==deltaIavg) {
......@@ -154,6 +176,9 @@ int main(int narg, char* arg[])
// puissance moyenne ds +/1 MHz @1420
fnorm=(1./(double)deltaIavg)/(double)(JFmax21-JFmin21+1);
for(int k=0; k<8; k++) xnt[2+8+k]*=fnorm;
// puissance moyenne en bande large B
fnorm=(1./(double)deltaIavg)/(double)(JFmaxB-JFminB+1);
for(int k=0; k<8; k++) xnt[2+24+k]*=fnorm;
nt.Fill(xnt);
I=0; cntnt++;
}
......
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