Skip to content

GitLab

Commit a906c88f authored by Magneville's avatar Magneville
Browse files

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

parents ccbb88ca 47dbe520
Hide whitespace changes
Inline Side-by-side
Showing with 475 additions and 21 deletions
anap4_exemple.d
View file @ a906c88f
......@@ -26,6 +26,9 @@ phcor PhaseCorrParamFileName
#### Time - frequency maps
## @tfm TFM_TimeBin TFM_FreqBin
@tfm 10 4
#### Right Ascension maps
## @ra RA_center_hours RA_width_hours RA_resol_minutes
@ra 12. 24. 1.
#### PrintLevel
# @prtlev PrintLevel PrintModulo
@prtlev 1 20
......
makefile
View file @ a906c88f
......@@ -10,13 +10,13 @@ MYINCLISTHERE = p4autils.h visip4reader.h visp4winreader.h p4gnugain.h p4gvcor.h
MYOLISTHERE = $(OBJ)/p4autils.o $(OBJ)/visip4reader.o $(OBJ)/visp4winreader.o $(OBJ)/p4gnugain.o $(OBJ)/p4gvcor.o $(OBJ)/p4phcor.o
# Define our target list
all : rdvisip4 visi2ntac visi2dtacx visi2tmfreq p4conv2fits msvis2dt visiavg filt_blind ckconvphases
all : rdvisip4 visi2ntac visi2dtacx visi2tmfreq p4conv2fits msvis2dt visiavg filt_blind ckconvphases vis2ra
clean :
rm -f $(EXE)/rdvisip4 $(EXE)/visi2ntac $(EXE)/visi2dtacx $(EXE)/visi2tmfreq $(EXE)/p4conv2fits $(EXE)/msvis2dt $(EXE)/visiavg \
$(EXE)/filt_blind $(EXE)/ckconvphases
$(EXE)/filt_blind $(EXE)/ckconvphases $(EXE)/vis2ra
rm -f $(OBJ)/rdvisip4.o $(OBJ)/visi2ntac.o $(OBJ)/visi2dtacx.o $(OBJ)/visi2tmfreq.o $(OBJ)/p4conv2fits.o $(OBJ)/msvis2dt.o $(OBJ)/visiavg.o \
$(OBJ)/filt_blind.o $(OBJ)/ckconvphases.o
$(OBJ)/filt_blind.o $(OBJ)/ckconvphases.o $(OBJ)/vis2ra.o
rm -f $(MYOLISTHERE)
depend :
......@@ -59,6 +59,17 @@ $(EXE)/visiavg : $(OBJ)/visiavg.o $(MYOLISTHERE)
$(OBJ)/visiavg.o : visiavg.cc $(MYINCLISTHERE)
$(CXXCOMPILE) -o $(OBJ)/visiavg.o visiavg.cc
## Calcul de cartes en ra (Right Ascension) a partir de donnees de visibilites
vis2ra : $(EXE)/vis2ra
echo '---vis2ra made'
$(EXE)/vis2ra : $(OBJ)/vis2ra.o $(MYOLISTHERE)
$(CXXLINK) -o $(EXE)/vis2ra $(OBJ)/vis2ra.o $(MYOLISTHERE) $(SOPHYAEXTSLBLIST)
$(OBJ)/vis2ra.o : vis2ra.cc $(MYINCLISTHERE)
$(CXXCOMPILE) -o $(OBJ)/vis2ra.o vis2ra.cc
## filtrage median dun canal
filt_blind : $(EXE)/filt_blind
echo '---filt_blind made'
......
p4autils.cc
View file @ a906c88f
......@@ -200,7 +200,8 @@ P4AnaParams::P4AnaParams()
fgserall_(false), fgtmsel_(false), inwsz_(1), fgreorderfreq_(false),
gain_gnu_file_(""), gain_var_file_(""), gain_blind_filt_file_(""), phase_corr_param_file_(""),
outfile_("p4a.ppf"), fitsoutfile_(""), prtlev_(0), prtmodulo_(1),
Nmean_(10), fgTFM_(false), TFMtimebin_(30), TFMfreqbin_(4), fgdtable_(false)
Nmean_(10), fgTFM_(false), TFMtimebin_(30), TFMfreqbin_(4), fgdtable_(false),
fgRA_(false), RAcenter_(12.), RAwidth_(24.), RAresol_mn_(1.)
{
}
......@@ -213,7 +214,9 @@ P4AnaParams::P4AnaParams(P4AnaParams const & a)
phase_corr_param_file_(a.phase_corr_param_file_),
outfile_(a.outfile_), fitsoutfile_(a.fitsoutfile_), prtlev_(a.prtlev_), prtmodulo_(a.prtmodulo_),
Nmean_(a.Nmean_), fgTFM_(a.fgTFM_), TFMtimebin_(a.TFMtimebin_), TFMfreqbin_(a.TFMfreqbin_),
fgdtable_(a.fgdtable_), fbands_(a.fbands_), lastargs_(a.lastargs_)
fgdtable_(a.fgdtable_), fbands_(a.fbands_),
fgRA_(a.fgRA_), RAcenter_(a.RAcenter_), RAwidth_(a.RAwidth_), RAresol_mn_(a.RAresol_mn_),
lastargs_(a.lastargs_)
{
}
......@@ -228,7 +231,9 @@ P4AnaParams& P4AnaParams::operator = (P4AnaParams const & a)
outfile_=a.outfile_; fitsoutfile_=a.fitsoutfile_;
prtlev_=a.prtlev_; prtmodulo_=a.prtmodulo_;
Nmean_=a.Nmean_; fgTFM_=a.fgTFM_; TFMfreqbin_=a.TFMfreqbin_; TFMtimebin_=a.TFMtimebin_;
fgdtable_=a.fgdtable_; fbands_=a.fbands_; lastargs_=a.lastargs_;
fgdtable_=a.fgdtable_; fbands_=a.fbands_;
fgRA_=a.fgRA_; RAcenter_=a.RAcenter_; RAwidth_=a.RAwidth_; RAresol_mn_=a.RAresol_mn_;
lastargs_=a.lastargs_;
return *this;
}
......@@ -308,6 +313,12 @@ int P4AnaParams::DecodeArgs(int narg, const char* arg[])
sscanf(arg[2],"%lg,%lg",&f0,&df); fbands_.push_back(P4FreqBand(f0,df)); fgdtable_=true;
arg+=2; narg-=2; lastargs_.clear();
}
else if (fbo=="-ra") {
if (narg<2) { cout << "P4AnaParams::DecodeArgs missing/bad argument, -h for help " << endl; return -1; }
int ja,jb;
sscanf(arg[2],"%lg,%lg,%lg",&RAcenter_,&RAwidth_,&RAresol_mn_); fgRA_=true;
arg+=2; narg-=2; lastargs_.clear();
}
else if (fbo=="-prt") {
if (narg<2) { cout << "P4AnaParams::DecodeArgs missing/bad argument, -h for help " << endl; return -1; }
sscanf(arg[2],"%d,%d",&prtlev_,&prtmodulo_);
......@@ -399,7 +410,11 @@ int P4AnaParams::ReadDatacardFile(const char * dcfilename)
cnt++;
}
if (dc.HasKey("tfm")) { // @tfm TFM_TimeBin TFM_FreqBin
TFMtimebin_=dc.IParam("tfm",0,1); TFMfreqbin_=dc.IParam("tfm",1,1);
TFMtimebin_=dc.IParam("tfm",0,1); TFMfreqbin_=dc.IParam("tfm",1,1); fgTFM_=true;
cnt++;
}
if (dc.HasKey("ra")) { // @ra RA_center_hours RA_width_hours RA_resol_minutes
RAcenter_=dc.DParam("ra",0,12.); RAwidth_=dc.DParam("ra",1,24.); RAresol_mn_=dc.DParam("ra",2,1.); fgRA_=true;
cnt++;
}
if (dc.HasKey("prtlev")) { // @prtlev PrintLevel PrintModulo
......@@ -418,7 +433,7 @@ int P4AnaParams::UsageOptions()
<< " [-dcf DataCardFileName ] [-in5 InputPath_BAO5] [-in6 InputPath_BAO6] \n"
<< " [-inrange Imin,Imax[,Istep] ] [-inall] [-inseltm YYYY-MM-DDThh:mm:ss duration_minutes] \n"
<< " [-o outfilename] [-gnu gain_filename] [-nmean N] [-reorderfreq] \n"
<< " [-tfm timebin,freqbin] [-freq f0,df] [-prt level] \n"
<< " [-tfm timebin,freqbin] [-ra RAcenterHours,RAwidthHours,RAresolMinutes] [-freq f0,df] [-prt level] \n"
<< " -dcf DatacardFileName : read and decode parameters from datacard file \n"
<< " -in5/6 InputPath_BAO5/6 : input file path for BAO5 and BAO6 \n"
<< " -inrange Imin,Imax[,Istep] : define range of visibilite matrices to be processed (Imin<=I<=Imax with step Istep) \n"
......@@ -435,6 +450,7 @@ int P4AnaParams::UsageOptions()
<< " -nmean N: compute average (and sigma...) using N visibility matrices \n"
<< " -reorderfreq: activate frequency reordering in VisiP4Reader, suitable for FFT firmware \n"
<< " -tfm timebin,freqbin : time-frequency maps, averaging timebin=nb of visi matrices, freqbin frequency components \n"
<< " -ra RAcenterHours,RAwidthHours,RAresolMinutes: define RA (right ascension) map center,total width (hours), resolution (minutes) \n"
<< " -freq f0,df : definition of a frequency band centered on f0, with total width df (in MHz) \n"
<< " several frequency bands can be specified by mutiple -freq options \n"
<< " -prt level,modulo: define print level and print modulo count \n"
......@@ -468,6 +484,8 @@ ostream& P4AnaParams::Print(ostream& os) const
os<<" Gain variation : parameters in "<< gain_var_file_<< "\n filtered vectors in "<< gain_blind_filt_file_ <<endl;
os<<" Phase correction : parameters in "<< phase_corr_param_file_ <<endl;
os<<" Nmean="<<Nmean_<<" TFM: TimeBin="<<TFMtimebin_<<" FreqBin="<<TFMfreqbin_<<endl;
os<<" Right Ascension: center_="<<RAcenter_<<" Width="<<RAwidth_<<" (hours) Resolution="<<RAresol_mn_
<<"("<<RAcenter_-RAwidth_/2<<"<ra<"<<RAcenter_-RAwidth_/2<<")"<<endl;
if (fbands_.size()==0) os<<" No frequency band defined ..."<<endl;
else for(size_t i=0; i<fbands_.size(); i++) os<<"FreqBand["<<i<<"]: "<<fbands_[i]<<endl;
return os;
......@@ -493,3 +511,43 @@ double P4Coords::RAFromTimeTU(TimeStamp const & tu){
return (ra);
}
//------------------------------------------------------
//----------------- class P4RAMapUtil ------------------
//------------------------------------------------------
P4RAMapUtil::P4RAMapUtil(P4AnaParams const& param)
: RAcenter_(param.RAcenter_), RAwidth_(param.RAwidth_), RAresol_mn_(param.RAresol_mn_)
{
if ((RAcenter_<0.)||(RAcenter_>=24.)||(RAwidth_<=1.e-19)||(RAwidth_>24.)||(RAresol_mn_<1.e-19)||(RAresol_mn_>360.))
throw ParmError("P4RAMapUtil::P4RAMapUtil(P4AnaParams const&) Bad RAcenter_ or RAwidth_ or RAresol_mn_ parameters");
invres_hours_ = 60./RAresol_mn_;
size_ = (sa_size_t)(RAwidth_*invres_hours_);
size_24hours_ = (sa_size_t)(24.*invres_hours_);
cout << " P4RAMapUtil: Resol(min)="<<RAresol_mn_<<" size="<< size_<<" size for 24 hours="<<size_24hours_<<endl;
}
P4RAMapUtil::P4RAMapUtil(double RAcenter, double RAwidth, double RAresol_mn)
: RAcenter_(RAcenter), RAwidth_(RAwidth), RAresol_mn_(RAresol_mn)
{
if ((RAcenter_<0.)||(RAcenter_>=24.)||(RAwidth_<=1.e-19)||(RAwidth_>24.)||(RAresol_mn_<1.e-19)||(RAresol_mn_>360.))
throw ParmError("P4RAMapUtil::P4RAMapUtil(RAcenter,...) Bad RAcenter_ or RAwidth_ or RAresol_mn_ parameters");
invres_hours_ = 60./RAresol_mn_;
size_ = (sa_size_t)(RAwidth_*invres_hours_);
size_24hours_ = (sa_size_t)(24.*invres_hours_);
cout << " P4RAMapUtil: Resol(min)="<<RAresol_mn_<<" size="<< size_<<" size for 24 hours="<<size_24hours_<<endl;
}
sa_size_t P4RAMapUtil::getMapIndex(double ra)
{
double delra = ra-RAcenter_;
sa_size_t idx = (sa_size_t)(delra*invres_hours_)+size_/2;
//DBG cout << " delra="<<delra<<" idx="<<idx<<endl;
if (idx<0) { // on ajoute 24 heures a ra et on recommece
idx += size_24hours_;
}
else if (idx>=size_) { // on soustrait 24 heures a ra et on retente
idx -= size_24hours_;
}
if ((idx>=0)&&(idx<size_)) return idx;
return -1;
}
p4autils.h
View file @ a906c88f
......@@ -206,6 +206,9 @@ public:
// definition des bandes de frequences pour les DataTables
bool fgdtable_;
vector<P4FreqBand> fbands_;
// Definition de la carte en ra (Right Ascension) - heures decimales, sauf pour la resolution, en minutes
bool fgRA_; // True , on a defini les limites de la carte en RA
double RAcenter_, RAwidth_, RAresol_mn_;
// les arguments en fin de ligne de commande
vector<string> lastargs_;
};
......@@ -227,4 +230,19 @@ public :
//! static method Converting a TU time into RA (Right Ascension sky coordinates) in decimal degrees
static double RAFromTimeTU(TimeStamp const & tu);
};
//! Class to facilitate projection into RA map
class P4RAMapUtil {
public:
P4RAMapUtil(P4AnaParams const& param);
P4RAMapUtil(double RAcenter, double RAwidth=24., double RAresol_mn=1.);
inline sa_size_t getMapSize() { return size_; }
//! return the index in map for ra (in decimal hours) - return negative if outside map
sa_size_t getMapIndex(double ra);
double RAcenter_, RAwidth_, RAresol_mn_;
double invres_hours_;
sa_size_t size_;
sa_size_t size_24hours_;
};
#endif
vis2ra.cc 0 → 100644
View file @ a906c88f
// Utilisation de SOPHYA pour faciliter les tests ...
#include "sopnamsp.h"
#include "machdefs.h"
/* ----------------------------------------------------------
Projet BAORadio/PAON4 - (C) LAL/IRFU
visiavg: programme de lecture des fichiers matrices de
visibilites de PAON4, calcul cartes de visibilites en
fonction de ra (Right Ascension)
O. Perdereau, R.Ansari - LAL , Janvier 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 "tvector.h"
#include "fioarr.h"
// #include "tarrinit.h"
#include "ntuple.h"
#include "datatable.h"
#include "histinit.h"
#include "matharr.h"
#include "timestamp.h"
#include <utilarr.h>
// include sophya mesure ressource CPU/memoire ...
#include "resusage.h"
#include "ctimer.h"
#include "timing.h"
// include lecteur de fichiers visibilites
#include "visp4winreader.h"
#include "fitsioserver.h"
#include "fiosinit.h"
int Usage(void);
int Usage(void)
{
cout << " --- vis2ra.cc : Reads PPF files produced by mfacq time-frequency \n"
<< " and produces maps of visibilies as a function of ra (Right Ascension) \n"
<< " Usage: vis2ra [-arguments] \n" << 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();
FitsIOServerInit();
P4AnaParams params;
params.DecodeArgs(narg, arg);
string outfile = params.outfile_;
if (outfile.length()<1) outfile = "visavg.ppf";
string fitsoutfile = params.fitsoutfile_;
if (fitsoutfile.length()>=1) {
fitsoutfile = "!"+fitsoutfile ; // adds '!' ?
}
int deltaIavg = params.TFMtimebin_;
sa_size_t TFMfbin = params.TFMfreqbin_;
int prtlev = params.prtlev_;
bool FgTFMAC = true;
bool FgTFMCX = true;
string desctfmap;
bool FgTFM = params.fgTFM_; // true -> create time-frequency maps
params.Print(cout);
cout <<"vis2ra/Info: Path BAO5:"<<params.inpath5_<<" BAO6:"<<params.inpath6_<<"\n"
<<"fgreorderfreq="<<params.fgreorderfreq_<<"\n"
<<" DeltaIAvg="<<deltaIavg<<"\n"
<<"outfile="<<outfile<<" PrtLev="<<prtlev<<endl;
P4RAMapUtil ram(params);
P4AVisiNumEncoder visiencod;
vector<sa_size_t> KVAC = visiencod.getAllAutoCor();
vector<sa_size_t> KVCXHH = visiencod.getAllHCrossCor();
cout << " List of AutoCorrelation rows:"<<endl;
for(size_t k=0; k<KVAC.size(); k++) {
cout << "KVAC["<<k<<"]="<<KVAC[k]<<" ->"<<visiencod.Convert2VisiName(KVAC[k])<<endl;
}
cout << " List of HH X-cor rows:"<<endl;
for(size_t k=0; k<KVCXHH.size(); k++) {
cout << "KVCXHH["<<k<<"]="<<KVCXHH[k]<<" ->"<<visiencod.Convert2VisiName(KVCXHH[k])<<endl;
}
// ---
HiStatsInitiator _inia;
int rc = 0;
try {
ResourceUsage resu;
VisiP4WindowReader wreader(params);
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;
bool fgok=true;
// vecteur de noms
vector <string> ext_names;
// un vecteur avec les temps
// TVector< double > timevec(wreader.getTotalNbWindows()/deltaIavg);
TMatrix< complex<r_4> > vismtx;
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;
//----- Count map for projection into (RA,Freq) maps
TArray< r_4 > racntmap;
//----- 6 H-H cross-cor TimeFrequency maps
vector< TArray< complex<r_4> > > vram;
//----- 6 H-H cross-cor TimeFrequency maps for the variances of real and imag parts
vector< TArray< r_4 > > vram_rp_sq;
vector< TArray< r_4 > > vram_ip_sq;
//----- 8 auto-corr TimeFrequency maps
vector< TArray< r_4 > > vramac;
//----- 8 auto-corr TimeFrequency variance maps
vector< TArray< r_4 > > vramac_sq;
sa_size_t ramSX, ramSY;
while (fgok) {
//reads next visibility matrix window
fgok = wreader.Shift();
if (!fgok) break;
vismtx = wreader.getAverageVisMtx(cfdate);
if (cnt==0) { //resizing matrices for sum of auto-correlations and sum of 6 cross-correlations
ramSX=ram.getMapSize();
ramSY=vismtx.NCols()/TFMfbin;
racntmap.SetSize2D(ramSX, ramSY, true); // true pour forcer la mise a zero
//allocating 8 Auto-Corr RA-frequency maps
cout<<"vis2ra/Info: allocating 8 AutoCor rightascension-Frequency maps : RA->NX="<<ramSX<<" x Freq->NY="<<ramSY<<endl;
for(int k=0; k<8; k++) vramac.push_back( TArray< r_4 >(ramSX, ramSY) );
cout <<" and 8 for the the variance maps "<<endl;
for(int k=0; k<8; k++) vramac_sq.push_back( TArray< r_4 >(ramSX, ramSY) );
//allocating 6 Cross-Corr H-H time-frequency maps
cout<<"vis2ra/Info: allocating H-H cross-cor RA-Frequency maps : RA->NX="<<ramSX<<" x Freq->NY="<<ramSY<<endl;
for(int k=0; k<6; k++) vram.push_back( TArray< complex<r_4> >(ramSX, ramSY) );
cout << "and the 2x6 for squares of rela & imag parts " << endl;
for(int k=0; k<6; k++) vram_rp_sq.push_back( TArray< r_4 >(ramSX, ramSY) );
for(int k=0; k<6; k++) vram_ip_sq.push_back( TArray< r_4 >(ramSX, ramSY) );
// recupere le jour de depart @ 0h
datestart = TimeStamp(cfdate.DaysPart(),0.);
}
dateobs=cfdate;
double racourant = P4Coords::RAFromTimeTU(dateobs);
sa_size_t idxra = ram.getMapIndex(racourant);
if (idxra < 0) continue;
// On est dans la carte en RA , on remplit d'abord la carte de count
for(sa_size_t jy=0; jy<ramSY; jy++) racntmap(idxra,jy) += 1.;
//---- On s'occupe d'abord des autocorrelations P1 ... P8
for(size_t k=0; k<KVAC.size(); k++) { // Loop over the 8 auto-correlations
TVector<r_4> vac = real(vismtx.Row(k));
TVector<r_4> vacsq = vac.SquareElt();
TArray< r_4 > & ramap = vramac[k];
TArray< r_4 > & ramap_sq = vramac_sq[k];
for(sa_size_t jy=0; jy<ramSY; jy++) { // frequency binning
ramap(idxra, jy) += vac( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
ramap_sq(idxra, jy) += vacsq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
}
} //----- end of loop over the 8 AutoCor
//---- On s'occupe des 6 cross-correlations 1H-2H ... 3H-4H
for(size_t k=0; k<KVCXHH.size(); k++) { // loop over the 6 Xcor
TVector< complex<r_4> > vcx = vismtx.Row(k);
TVector<r_4> vcxprsq = (real(vcx)).SquareElt();
TVector<r_4> vcxpisq = (imag(vcx)).SquareElt();
TArray< complex<r_4> > & ramap = vram[k];
TArray< r_4 > & ramapsqpr = vram_rp_sq[k];
TArray< r_4 > & ramapsqpi = vram_ip_sq[k];
for(sa_size_t jy=0; jy<ramSY; jy++) {
ramap(idxra, jy) = vcx( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
ramapsqpr(idxra, jy) = vcxprsq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
ramapsqpi(idxra, jy) = vcxpisq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
}
} //----- end of loop over the 6 Xcor
// double tdif = cfdate.TimeDifferenceSeconds(cfdate,dateobs)/2.;
// timevec(TFMtmidx) = dateobs.TimeDifferenceSeconds(dateobs.ShiftSeconds (tdif ),datestart); // centre du bin
// ... done
cnt++;
if ( (cnt>0) && (cnt%50==0) ) {
// cout << " vis2ra/Info: RA-map fill, VisMtxCount="<<cnt<< " /Max="<<Imax<<" DateObs="<<dateobs<<endl;
}
} // Fin de la boucle de lecture
cout<<"vis2ra/Info: count="<<cnt<<" visimtx read "<<endl;
// --- Sauvegarde cartes temps-frequence en fits
//FitsABTWriter * fbtw = NULL;
FitsInOutFile * fos = NULL ;
cout<<"vis2ra/Info: Normalising RA-maps ... " << endl;
for(size_t k=0; k<KVAC.size(); k++) { // loop over 8 Autocor
TArray< r_4 > & ramap = vramac[k];
TArray< r_4 > & ramap_sq = vramac_sq[k];
for(sa_size_t ix=0; ix<ramSX; ix++) {
for(sa_size_t jy=0; jy<ramSY; jy++) {
if (racntmap(ix,jy) < 1.e-9) continue;
r_4 fnorm = 1./racntmap(ix,jy);
ramap(ix, jy) *= fnorm;
ramap_sq(ix, jy) *= fnorm;
}
}
} // end of loop over 8 Autocor
for(size_t k=0; k<KVCXHH.size(); k++) { // loop over the 6 Xcor
TArray< complex<r_4> > & ramap = vram[k];
TArray< r_4 > & ramapsqpr = vram_rp_sq[k];
TArray< r_4 > & ramapsqpi = vram_ip_sq[k];
for(sa_size_t ix=0; ix<ramSX; ix++) {
for(sa_size_t jy=0; jy<ramSY; jy++) {
if (racntmap(ix,jy) < 1.e-9) continue;
r_4 fnorm = 1./racntmap(ix,jy);
complex<r_4> znorm = complex<r_4>(fnorm,0.);
ramap(ix, jy) *= znorm;
ramapsqpr(ix, jy) *= fnorm;
ramapsqpi(ix, jy) *= fnorm;
}
}
} // End of loop over the 6 Xcor
if (fitsoutfile.length()>=1){
cout << " fitsoutfile" <<fitsoutfile<< endl;
fos = new FitsInOutFile(fitsoutfile, FitsInOutFile::Fits_Create);
}
POutPersist poram(outfile);
char bufnam[10];
int numkey =0;
// --- saving AutoCorr time-frequency maps
cout<<" vis2ra/Info: Saving 8 AutoCorr time-frequency maps to PPF file "<<outfile<<endl;
const char* ram_names[8]={"RAM_1H", "RAM_2H", "RAM_3H", "RAM_4H", "RAM_1V", "RAM_2V", "RAM_3V", "RAM_4V"};
const char* ramsq_names[8]={"VARRAM_1H", "VARRAM_2H", "VARRAM_3H", "VARRAM_4H", "VARRAM_1V", "VARRAM_2V", "VARRAM_3V", "VARRAM_4V"};
for(int k=0; k<8; k++) { // loop over the 8 AutoCorr
TArray< r_4 > & ramap = vramac[k];
poram << PPFNameTag(ram_names[k]) << ramap;
if (fos != NULL) {
ext_names.push_back(ram_names[k]);
(*fos)<< ramap;
}
TArray< r_4 > & ramapsq = vramac_sq[k];
ramapsq = ramapsq - ramap.SquareElt();
poram << PPFNameTag(ramsq_names[k]) << ramapsq;
if (fos != NULL) {
ext_names.push_back(ramsq_names[k]);
(*fos)<< ramapsq;
}
}
// --- renormalizing and saving H-H Cross-Corr time-frequency maps
cout<<" vis2ra/Info: Saving 6 H-H cross-corr time-frequency maps to PPF file "<<outfile<<endl;
const char* ramCC_names[6]={"RAM_1H2H", "RAM_1H3H", "RAM_1H4H", "RAM_2H3H", "RAM_2H4H", "RAM_3H4H"};
const char* vrramCC_names[6]={"RVARRAM_1H2H", "RVARRAM_1H3H", "RVARRAM_1H4H", "RVARRAM_2H3H", "RVARRAM_2H4H", "RVARRAM_3H4H"};
const char* viramCC_names[6]={"IVARRAM_1H2H", "IVARRAM_1H3H", "IVARRAM_1H4H", "IVARRAM_2H3H", "IVARRAM_2H4H", "IVARRAM_3H4H"};
for(int k=0; k<6; k++) { // loop over the 6 Xcor
TArray< complex<r_4> > & ramap = vram[k];
TArray< r_4 > & ramap_sqpr = vram_rp_sq[k];
TArray< r_4 > & ramap_sqpi = vram_ip_sq[k];
poram << PPFNameTag(ramCC_names[k]) << ramap;
if (fos != NULL) {
ext_names.push_back(string(ramCC_names[k])+"_real");
(*fos)<< real(ramap);
ext_names.push_back(string(ramCC_names[k])+"_imag");
(*fos)<< imag(ramap);
}
TArray< r_4 > ramapr = real(ramap);
TArray< r_4 > ramapi = imag(ramap);
ramap_sqpr -= ramapr.SquareElt() ;
ramapi = ramapi.MulElt(ramapi,ramapi) ;
ramap_sqpi -= ramapi.SquareElt() ;
poram << PPFNameTag(vrramCC_names[k]) << ramap_sqpr;
poram << PPFNameTag(viramCC_names[k]) << ramap_sqpi;
if (fos != NULL) {
ext_names.push_back(vrramCC_names[k]);
(*fos)<< ramap_sqpr;
ext_names.push_back(viramCC_names[k]);
(*fos)<< ramap_sqpi;
}
}
// --- FIN sauvegarde cartes temps-frequence
// resu.Update();
P4FreqBand myp4fre;
TVector <double> avg_freqs( myp4fre.getP4NbFreqChannels()/TFMfbin);
double frbase = myp4fre.freqstart_ + myp4fre.getP4FreqResolution()/2. ;
for (int kf=0 ; kf< myp4fre.getP4NbFreqChannels()/TFMfbin ; kf++,frbase += myp4fre.getP4FreqResolution()*TFMfbin )
avg_freqs(kf) = frbase ;
poram << PPFNameTag("FreqVec") << avg_freqs;
if (fos != NULL) {
ext_names.push_back("Frequences");
(*fos)<< avg_freqs ;
cout << " number of objs in fits "<< ext_names.size() << endl;
cout << ext_names << endl;
delete(fos);
}
cout << resu; // Update est fait lors du print
}
catch (PException& exc) {
cerr << " vis2ra.cc catched PException " << exc.Msg() << endl;
rc = 77;
}
catch (std::exception& sex) {
cerr << "\n vis2ra.cc std::exception :"
<< (string)typeid(sex).name() << "\n msg= "
<< sex.what() << endl;
rc = 78;
}
catch (...) {
cerr << " vis2ra.cc catched unknown (...) exception " << endl;
rc = 79;
}
cout << ">>>> vis2ra.cc ------- END ----------- RC=" << rc << endl;
return rc;
}
visp4winreader.cc
View file @ a906c88f
......@@ -9,24 +9,25 @@
/* --Methode-- */
VisiP4WindowReader::VisiP4WindowReader(vector<string> const& vpath, size_t wsz)
: wsz_(wsz) , p4vreaderp_(NULL), vtu_(wsz),
p4gnu_(NULL), p4gvf_(NULL), p4phc_(NULL),
fgown_p4gnu_(false), fgown_p4gvf_(false), fgown_p4phc_(false)
: wsz_(wsz) , p4vreaderp_(NULL), vtu_(wsz),
p4gnu_(NULL), p4gvf_(NULL), p4phc_(NULL),
fgown_p4gnu_(false), fgown_p4gvf_(false), fgown_p4phc_(false), fggpcor_(wsz)
{
cout << " VisiP4WindowReader - Creating Reader with window size= " << wsz_ << endl;
BaseArray::SetDefaultMemoryMapping(BaseArray::CMemoryMapping);
for(size_t ii=0; ii<wsz_; ii++) fggpcor_[ii]=false; cntfgcor_=0;
p4vreaderp_ = VisiP4ReaderBase::getReader(vpath);
}
/* --Methode-- */
VisiP4WindowReader::VisiP4WindowReader(P4AnaParams & params)
: wsz_(params.inwsz_) , p4vreaderp_(NULL), vtu_(params.inwsz_),
p4gnu_(NULL), p4gvf_(NULL), p4phc_(NULL),
fgown_p4gnu_(false), fgown_p4gvf_(false), fgown_p4phc_(false)
: wsz_(params.inwsz_) , p4vreaderp_(NULL), vtu_(params.inwsz_),
p4gnu_(NULL), p4gvf_(NULL), p4phc_(NULL),
fgown_p4gnu_(false), fgown_p4gvf_(false), fgown_p4phc_(false), fggpcor_(params.inwsz_)
{
cout << " VisiP4WindowReader(P4AnaParams & params) - Creating Reader with window size= " << wsz_ << endl;
BaseArray::SetDefaultMemoryMapping(BaseArray::CMemoryMapping);
for(size_t ii=0; ii<wsz_; ii++) fggpcor_[ii]=false; cntfgcor_=0;
vector<string> paths;
paths.push_back(params.inpath5_);
paths.push_back(params.inpath6_);
......@@ -105,6 +106,7 @@ void VisiP4WindowReader::setPhaseCorrector(P4PhaseCor& phc)
/* --Methode-- */
bool VisiP4WindowReader::Shift()
{
for(size_t ii=0; ii<wsz_; ii++) fggpcor_[ii]=false; cntfgcor_=0;
size_t k0 = 0;
TimeStamp dateobs;
double mttag;
......@@ -134,22 +136,26 @@ TMatrix< complex<r_4> > VisiP4WindowReader::getVisMtx(size_t k, bool fgcor)
if (k >= wsz_) throw RangeCheckError("VisiP4WindowReader::getVisMtx(k) k >= WindowSize");
TMatrix< complex<r_4> > vismtx = visarr_(Range::all(), Range::all(), Range((sa_size_t)k));
if (!fgcor || (!p4gnu_ && !p4gvf_ && !p4phc_)) return vismtx;
TMatrix< complex<r_4> > rvmx(vismtx, false); // copy the visibility matrix to avoid changing the array