Commit 2110933e authored by Reza  ANSARI's avatar Reza ANSARI
Browse files

1/ Adaptation des parametres (-inth) pour la lecture des voies thermometre et RFI

2/ Adaptation de la classe VisiP4WindowReader pour la lecture des voies Thermometre et RFI
3/ Ajout programme rdthermrfi.cc de lecture voies Thermometre et RFI et adaptation Makefile
                                          Reza , 18/07/2018
parent 5b1c1e2d
// Utilisation de SOPHYA pour faciliter les tests ...
#include "sopnamsp.h"
#include "machdefs.h"
/* ----------------------------------------------------------
Projet BAORadio/PAON4 - (C) LAL/IRFU 2008-2018
Programme de lecture des fichiers matrices de visibilites de PAON4
produits par mfacq pour les voies Thermometres et RFI (9,10)
R. Ansari Univ. Paris Sud, & LAL CNRS/IN2P3
Juillet 2018
---------------------------------------------------------- */
// include standard c/c++
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <iostream>
#include <string>
#include "pexceptions.h"
#include "array.h"
#include "histats.h"
#include "fitsioserver.h"
#include "resusage.h"
// include lecteur de fichiers visibilites
#include "p4autils.h"
#include "visp4winreader.h"
int Usage(void)
{
cout<<"--- rdthermrfi : Read PPF files for thermometer & RFI signal channnels and produce TF-Maps \n"<<endl;
cout<<" Usage: rdthermrfi [-options] \n";
cout<<" Specify -inth PathForThermometer or in datacards " <<endl;
P4AnaParams::UsageOptions();
cout<<endl;
return 1;
}
//----------------------------------------------------
int main(int narg, const char* arg[])
{
// --- Decoding parameters
if ((narg<2)||((narg>1)&&(strcmp(arg[1],"-h")==0))) return Usage();
P4AnaParams params;
params.DecodeArgs(narg, arg);
if (params.inpath_thermRFI_.length()<1) {
cout<<"rdthermrfi/ERROR no path specified for thermometer/RFI signals ... (rdthermrfi -h for usage)" <<endl;
return 2;
}
string outfile = params.outfile_;
if (outfile.length()<1) outfile = "thrfiavg.ppf";
string fitsoutfile = params.fitsoutfile_;
if (fitsoutfile.length()>=1) {
fitsoutfile = "!"+fitsoutfile ; // adds '!' ?
}
int deltaIavg = params.TFMtimebin_;
sa_size_t TFMfbin = params.TFMfreqbin_;
int rc = 0;
try {
FitsIOServerInit();
HiStatsInitiator _inia;
ResourceUsage resu;
VisiP4WindowReader wreader(params, true); // true -> Read Therometer/RFI
long Imin = wreader.getReader().getSerialFirst();
long Imax = wreader.getReader().getSerialLast();
long Istep = wreader.getReader().getSerialStep();
cout << "vis2ra/Info: processing visibility matrix serial/sequence number range "
<<Imin<<" <= seq <= " << Imax << " with step="<<Istep<<endl;
cout << " WindowSize="<<wreader.getWindowSize()<<" -> TotalNbWindows="<<wreader.getTotalNbWindows()<<endl;
//vismtx: visib matrice, meanmtx: mean FFT coef, vismtx_center: visib matrix corrigee de l'offset (normalisee avec NPAQSUM)
// un vecteur avec les temps
TVector< double > timevec(wreader.getTotalNbWindows()/deltaIavg);
TVector< double > ravec(wreader.getTotalNbWindows()/deltaIavg);
TMatrix< complex<r_4> > meanvismtx;
TMatrix< complex<r_4> > visum;
TimeStamp dateobs, cfdate,datestart;
TimeStamp dateorg(2015,1,1,12,0,0.); // Date origine 1 jan 2015
double mttag;
int cnt=0, cntnt=0, pcntnt=0;
int I=0;
//----- 2 auto-corr (Thermometer/RFI) TimeFrequency maps
TArray<r_4> TFM_therm;
TArray<r_4> TFM_RFI;
//---- 1 cross-correlation TFM (Thermometer <> RFI)
TArray< complex<r_4> > TFM_cross;
//---- for the time-freqency map filling
sa_size_t TFMtmidx=0;
sa_size_t tfmSX, tfmSY;
bool fgok=true;
while (fgok) {
//reads next visibility matrix window
fgok = wreader.Shift();
if (!fgok) break;
TMatrix< complex<r_4> > vismtx = wreader.getAverageVisMtx(cfdate);
if (cnt==0) { //resizing matrices for sum of auto-correlations and sum of 6 cross-correlations
visum = vismtx; meanvismtx=vismtx;
tfmSX=wreader.getTotalNbWindows()/deltaIavg;
tfmSY=vismtx.NCols()/TFMfbin;
TFM_therm.SetSize2D(tfmSX, tfmSY);
TFM_RFI.SetSize2D(tfmSX, tfmSY);
TFM_cross.SetSize2D(tfmSX, tfmSY);
dateobs=cfdate;
continue;
}
if (I==0) dateobs=cfdate; // start filling a new time bin
// Accumulating
visum += vismtx;
meanvismtx += vismtx;
I++; // incrementing DeltaTime counter
// we check that our time index did not go beyond the allocated array size (might not be necessary)
if ((I==deltaIavg)&&(TFMtmidx>=tfmSX)) { // Cela ne devrait pas arriver en principe
TFMtmidx++; I=0;
cout << "rdthermrfi/Warning: something wrong in the logic , (TFMtmidx="<<TFMtmidx<<") >= (tfmSX="<<tfmSX<<")"
<< " for read count="<<cnt<<endl;
}
else if (I==deltaIavg) { // Filling TFM maps
TVector<r_4> vac0 = real(visum.Row(0));
TVector<r_4> vac1 = real(visum.Row(1));
TVector< complex<r_4> > vacX = visum.Row(2);
for(sa_size_t jy=0; jy<tfmSY; jy++) { // frequency binning
TFM_therm(TFMtmidx, jy) = vac0( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
TFM_RFI(TFMtmidx, jy) = vac1( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
TFM_cross(TFMtmidx, jy) = vacX( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
}
double tdif = cfdate.TimeDifferenceSeconds(cfdate,dateobs)/2.;
timevec(TFMtmidx) = dateobs.TimeDifferenceSeconds(dateobs.ShiftSeconds (tdif ),datestart); // centre du bin
ravec(TFMtmidx) = P4Coords::RAFromTimeTU(dateobs.ShiftSeconds (tdif ));
TFMtmidx++;
// ... done
I=0; cntnt++;
visum = complex<r_4>(0.,0.);
}
cnt++;
if ((cnt>0)&&(cntnt%10==0)&&(cntnt>pcntnt)) {
cout<<"rdthermrfi/Info: TFM-Map fill cnt="<<cntnt<<" VisMtxCount="<<cnt
<<" /Max="<<wreader.getTotalNbWindows()<<" DateObs="<<dateobs<<endl;
pcntnt=cntnt;
}
} // End of redaing visibility matrices loop
cout<<"rdthermrfi/Info: count="<<cnt*wreader.getWindowSize()<<" Visibility Matrices read "<<endl;
if (cnt>0) {
meanvismtx /= complex<r_4>((r_4)cnt, (r_4)0.);
cout<<" rdthermrfi/Info: Saving vismtx_mean to PPF file "<<outfile<<endl;
POutPersist po(outfile);
po<<PPFNameTag("meanvis_Th_RFI")<<meanvismtx;
cout<<" rdthermrfi/Info: TFM_therm , TFM_RFI, TFM_cross (TFM_ThXRFI) to PPF file "<<outfile<<endl;
po<<PPFNameTag("TFM_therm")<<TFM_therm;
po<<PPFNameTag("TFM_RFI")<<TFM_RFI;
po<<PPFNameTag("TFM_ThXRFI")<<TFM_cross;
}
// resu.Update();
cout << resu; // Update est fait lors du print
}
catch (PException& exc) {
cerr << " rdthermrfi.cc catched PException " << exc.Msg() << endl;
rc = 77;
}
catch (std::exception& sex) {
cerr << "\n rdthermrfi.cc std::exception :"
<< (string)typeid(sex).name() << "\n msg= "
<< sex.what() << endl;
rc = 78;
}
catch (...) {
cerr << " rdthermrfi.cc catched unknown (...) exception " << endl;
rc = 79;
}
cout << ">>>> rdthermrfi.cc ------- END ----------- RC=" << rc << endl;
return rc;
}
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