Commit 12c40754 authored by Reza Ansari's avatar Reza Ansari
Browse files

MAJ programme visi2ntac : remplissage d'un datatable avec auto-correlations et...

MAJ programme visi2ntac : remplissage d'un datatable avec auto-correlations et cross-correlations pour 2 frequences - pas de binning ni en temps, ni en frequence, Reza 6/11/2017
parent 2564a2b4
......@@ -42,76 +42,52 @@
#include "p4gnugain.h"
//--------------------------- Fonctions de ce fichier -------------------
int Usage(bool fgshort=true);
// int DecodeArgs(int narg, char* arg[]);
/* --Fonction-- */
int Usage(bool fgea)
int Usage(void);
int Usage(void)
{
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 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"
<< " 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) {
cout << " rdvisip4 -h for detailed instructions" << endl;
return 1;
}
*/
cout << " --- visi2ntac.cc : Read PPF files produced by mfacq and fills an NTuple \n" << endl;
cout << " Usage: visi2ntac [-arguments] [-cg] [freq1=1393] [freq2=1420.4] \n"
<< " -cg : apply gain corrections \n"<<endl;
P4AnaParams::UsageOptions();
cout<< endl;
return 1;
}
//----------------------------------------------------
//----------------------------------------------------
int main(int narg, char* arg[])
int main(int narg, const char* arg[])
{
if ((narg>1)&&(strcmp(arg[1],"-h")==0)) return Usage(false);
int aoff = 0;
bool fgreorderfreq=false;
if ((narg>1)&&(strcmp(arg[1],"-reorderfreq")==0)) {
fgreorderfreq=true;
aoff = 1;
// --- Decoding parameters
if( (narg<2) || ((narg>1)&&(strcmp(arg[1],"-h")==0) ) ) return Usage();
P4AnaParams params;
params.DecodeArgs(narg, arg);
string outfile = params.outfile_;
if (outfile.length()<1) outfile = "visdt.ppf";
int Imin = params.Imin_, Imax = params.Imax_, Istep = params.Istep_;
bool fgcorgain = false;
double freq1 = 1393.;
double freq2 = 1420.4;
if (params.lastargs_.size()>0) {
size_t off=0;
if (params.lastargs_[0] == "-cg") { fgcorgain=true; off=1; }
if (params.lastargs_.size() > off) freq1 = atof(params.lastargs_[off].c_str());
if (params.lastargs_.size() > off+1) freq2 = atof(params.lastargs_[off+1].c_str());
}
if (narg<9+aoff) return Usage(true);
string path5 = arg[1+aoff];
string path6 = arg[2+aoff];
int Imin,Imax,Istep=1;
sscanf(arg[3+aoff],"%d,%d",&Imin,&Imax);
Istep=1;
string gainfile=arg[4+aoff];
string outfile=arg[5+aoff];
// frequency range definition
sa_size_t JFmin=1300,JFmax=1460;
int jf1,jf2;
if ((narg>6)&&(strcmp(arg[6+aoff],"-")!=0)) sscanf(arg[6+aoff],"%d,%d",&jf1,&jf2);
JFmin=jf1; JFmax=jf2;
sa_size_t JFminB=2130,JFmaxB=2290;
if ((narg>7)&&(strcmp(arg[7+aoff],"-")!=0)) sscanf(arg[7+aoff],"%d,%d",&jf1,&jf2);
JFminB=jf1; JFmaxB=jf2;
// time averaging interval definition, in term of visibility files
int deltaIavg=10;
if (narg>8) deltaIavg=atoi(arg[8+aoff]);
if (deltaIavg<1) deltaIavg=1;
int prtlev=1;
if (narg>9) prtlev=atoi(arg[9+aoff]);
cout << " visi2ntac/Info: Path BAO5:"<<path5<<" BAO6:"<<path6<<"\n"
<< " Imin,max,step="<<Imin<<","<<Imax<<","<<Istep<<" OutFile:"<<outfile<<"\n"
<< JFmin<<" <=NumFreq<= "<<JFmax<<" "<<JFminB<<" <=NumFreqB<= "<<JFmaxB
<<" DeltaIAvg="<<deltaIavg<<" PrtLev="<<prtlev<<endl;
P4gnuGain p4g(gainfile);
Matrix const & gains = p4g.getGainMatrix();
// 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
cout << " pjHI +/-1 MHz @1420 MHz " << JFmin21<<" <=NumFreq<= "<<JFmax21
<< " pjHI +/-5 MHz @1420 MHz " << JFmin21_5<<" <=NumFreq<= "<<JFmax21_5 << endl;
P4FreqBand fb1(freq1,0.);
sa_size_t JF1 = fb1.getFirstFreqChannel();
P4FreqBand fb2(freq1,0.);
sa_size_t JF2 = fb1.getFirstFreqChannel();
params.Print(cout);
cout <<"visi2ntac/Info: Path BAO5:"<<params.inpath5_<<" BAO6:"<<params.inpath6_<<"\n"
<<"fgreorderfreq="<<params.fgreorderfreq_<<"\n"
<<"Imin,max,step="<<Imin<<","<<Imax<<","<<Istep
<<"outfile="<<outfile<<endl;
cout << (fgcorgain ? " APPLY GainCorrection " : " NO GainCorr") << " Freq1="<<freq1<<" MHz -> JF1="<<JF1
<<" Freq1="<<freq1<<" MHz -> JF1="<<JF1<<endl;
P4gnuGain * p4gp=NULL;
if (fgcorgain) p4gp = new P4gnuGain(params.gain_gnu_file_);
HiStatsInitiator _inia;
// TArrayInitiator _inia;
......@@ -119,90 +95,118 @@ int main(int narg, char* arg[])
int rc = 0;
try {
ResourceUsage resu;
const char* dtnames[51]={ "k","day","time",
"ac1h","ac2h","ac3h","ac4h","ac1v","ac2v","ac3v","ac4v",
"ac1hB","ac2hB","ac3hB","ac4hB","ac1vB","ac2vB","ac3vB","ac4vB",
"cx1h2h","cx1h3h","cx1h4h","cx2h3h","cx2h4h","cx3h4h",
"cx1h2hB","cx1h3hB","cx1h4hB","cx2h3hB","cx2h4hB","cx3h4hB",
"cx1v2v","cx1v3v","cx1v4v","cx2v3v","cx2v4v","cx3v4v",
"cx1v2vB","cx1v3vB","cx1v4vB","cx2v3vB","cx2v4vB","cx3v4vB",
"cx1h4v","cx2h4v","cx3h4v","cx4h4v",
"cx1h4vB","cx2h4vB","cx3h4vB","cx4h4vB" };
DataTable dt(384);
dt.AddIntegerColumn(dtnames[0]);
dt.AddIntegerColumn(dtnames[1]);
dt.AddFloatColumn(dtnames[2]);
for(int kk=0; kk<8; kk++)
dt.AddFloatColumn(dtnames[3+kk]);
for(int kk=0; kk<24; kk++)
dt.AddComplexColumn(dtnames[3+8+kk]);
for(int kk=0; kk<8; kk++)
dt.AddComplexColumn(dtnames[3+8+24+kk]);
P4AVisiNumEncoder p4venc;
std::vector<sa_size_t> kac = p4venc.getAllAutoCor();
std::vector<sa_size_t> kcxh = p4venc.getAllHCrossCor();
std::vector<sa_size_t> kcxv = p4venc.getAllVCrossCor();
std::vector<sa_size_t> kcxh4v(4);
kcxh4v[0] = p4venc.Convert2VisiNumber("1H4V");
kcxh4v[1] = p4venc.Convert2VisiNumber("2H4V");
kcxh4v[2] = p4venc.Convert2VisiNumber("3H4V");
kcxh4v[3] = p4venc.Convert2VisiNumber("4H4V");
vector<string> paths;
paths.push_back(path5);
paths.push_back(path6);
const char* ntnames[35]={"k","day","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",
"p1B","p2B","p3B","p4B","p5B","p6B","p7B","p8B" };
NTuple nt(35,ntnames,384,false); // float numbers have enough precision for us
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);
paths.push_back(params.inpath5_);
paths.push_back(params.inpath6_);
VisiP4ReaderBase * reader = VisiP4ReaderBase::getReader(paths);
VisiP4ReaderBase & vreader = (*reader);
vreader.setFreqReordering(fgreorderfreq);
vreader.setFreqReordering(params.fgreorderfreq_);
if (!params.fgserall_ && !params.fgtmsel_) {
cout << " vreader.SelectSerialNum(Imin="<<Imin<<" ,Imax="<<Imax<<" ,Istep="<<Istep<<")"<<endl;
vreader.SelectSerialNum(Imin,Imax,Istep);
vreader.setPrintLevel(prtlev);
}
else if (params.fgserall_) {
cout << " vreader.SelectAll() ... " << endl;
vreader.SelectAll();
}
else {
TimeStamp tustart = params.tmsel_tu_; tustart.ShiftSeconds(-params.tmsel_duration_*30.);
TimeStamp tuend = params.tmsel_tu_; tuend.ShiftSeconds(params.tmsel_duration_*30.);
cout << " vreader.SelectTimeFrame(TUStart="<<tustart.ToString()<<" ,TUEnd="<<tuend.ToString()
<<" ,Istep="<<Istep<<")"<<endl;
vreader.SelectTimeFrame(tustart, tuend, Istep);
}
vreader.setPrintLevel(params.prtlev_);
Imin = vreader.getSerialFirst(); Imax = vreader.getSerialLast(); Istep = vreader.getSerialStep();
cout << "visi2dtacx/Info: processing visibility matrix serial/sequence number range "
<<Imin<<" <= seq <= " << Imax << " with step="<<Istep<<endl;
bool fgok=true;
TMatrix< complex<r_4> > vismtx;
TMatrix< complex<r_4> > acsum;
TimeStamp dateobs;
TimeStamp dateorg(2015,1,1,12,0,0.); // Date origine 1 jan 2015
TimeStamp dateorg(2017,9,1,12,0,0.); // Date origine 1 sep 2017
int cnt=0;
// Getting en empty row_ptr
DataTableRowPtr rowp = dt.EmptyRowPtr();
TMatrix< complex<r_4> > vismtx;
double mttag;
int cnt=0, cntnt=0, pcntnt=0;
int I=0;
while (fgok) {
fgok=vreader.ReadNext(vismtx, dateobs, mttag);
if (fgok) {
if (cnt==0) acsum.SetSize(8, vismtx.NCols());
if (I==0) {
xnt[0]=cnt; xnt[1]=dateobs.DaysPart()-dateorg.DaysPart(); xnt[2]=dateobs.SecondsPart();
if (prtlev>0)
cout << "visi2ntac/Info: dateobs="<<dateobs<<" SecondsPart()="<<dateobs.SecondsPart()<<endl;
for(int k=0; k<32; k++) xnt[3+k]=0.;
acsum = complex<r_4>(0.,0.);
if (! fgok) break;
if (fgcorgain) p4gp->applyGain(vismtx);
// remplissage datatable
dt.NextRowPtr(rowp);
rowp(0)=(int_4)cnt;
rowp(1)=(int_4)(dateobs.DaysPart()-dateorg.DaysPart());
rowp(2)=(r_4)dateobs.SecondsPart();
size_t off = 3;
for(size_t k=0; k<kac.size(); k++) {
rowp(off+k) = (float)(vismtx(kac[k], JF1).real());
rowp(off+k+kac.size()) = (float)(vismtx(kac[k], JF2).real());
}
for(int k=0; k<8; k++) acsum.Row(k) += vismtx.Row(KVAC[k]);
I++;
if (I==deltaIavg) {
Vector vac(acsum.NCols());
for(int k=0; k<8; k++) {
// correct for gain(nu)
for(sa_size_t c=0; c<vac.Size(); c++) vac(c)=acsum(k,c).real()*gains(k,c);
xnt[3+k] += vac(freqs).Sum(); // integration en bande large
xnt[3+8+k] += vac(freqs21).Sum(); // integration a +/- 1 MHz @ 1420 MHz
xnt[3+16+k] += vac(freqs21_5).Sum(); // integration a +/- 5 MHz @ 1420 MHz
xnt[3+24+k] += vac(freqsB).Sum(); // integration en bande large B
off += (2*kac.size());
for(size_t k=0; k<kcxh.size(); k++) {
rowp(off+k) = complex<float>(vismtx(kcxh[k], JF1));
rowp(off+k+kcxh.size()) = complex<float>(vismtx(kcxh[k], JF2));
}
// puissance moyenne dans la bande large
double fnorm=(1./(double)deltaIavg)/(double)(JFmax-JFmin+1);
for(int k=0; k<8; k++) xnt[3+k]*=fnorm;
// puissance moyenne ds +/5 MHz @1420 , exclu +/- 1 MHz centre a 1420
fnorm=(1./(double)deltaIavg)/(double)((JFmax21_5-JFmin21_5)-(JFmax21-JFmin21));
for(int k=0; k<8; k++) xnt[3+16+k] = (xnt[3+16+k]-xnt[3+8+k])*fnorm;
// puissance moyenne ds +/1 MHz @1420
fnorm=(1./(double)deltaIavg)/(double)(JFmax21-JFmin21+1);
for(int k=0; k<8; k++) xnt[3+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[3+24+k]*=fnorm;
nt.Fill(xnt);
I=0; cntnt++;
}
cnt++;
if ((cntnt>0)&&(cntnt%10==0)&&(cntnt>pcntnt)) {
cout << " visi2ntac/Info: fill-count="<<cntnt<<" DateObs="<<dateobs<<endl;
pcntnt=cntnt;
off += (2*kcxh.size());
for(size_t k=0; k<kcxv.size(); k++) {
rowp(off+k) = complex<float>(vismtx(kcxv[k], JF1));
rowp(off+k+kcxv.size()) = complex<float>(vismtx(kcxv[k], JF2));
}
off += (2*kcxv.size());
for(size_t k=0; k<kcxh4v.size(); k++) {
rowp(off+k) = complex<float>(vismtx(kcxh4v[k], JF1));
rowp(off+k+kcxh4v.size()) = complex<float>(vismtx(kcxh4v[k], JF2));
}
cnt++;
if (cnt%50==0)
cout << " visi2ntac/Info: fill-count="<<cnt<<" DateObs="<<dateobs<<" /Max="<<(Imax-Imin+1)/Istep<<endl;
}
cout << " visi2ntac/Info: count="<<cnt<<" visimtx read "<<endl;
cout << nt;
cout<<" rdvisip4/Info: Saving auto-corr=f(time) ntuple to PPF file "<<outfile<<endl;
cout << dt;
cout<<" rdvisip4/Info: Saving dt DataTable to PPF file "<<outfile<<endl;
POutPersist po(outfile);
po<<nt;
po<<dt;
// resu.Update();
cout << resu; // Update est fait lors du print
}
......
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