Commit 3a18bfad authored by Reza  ANSARI's avatar Reza ANSARI
Browse files

1/ Modification de la classe P4gvCor afin de prendre en charge la correction du gain global

(normalisation)
2/ optimisation de l'application du gain G(t) - on peut instancier la classe avec un
constructeur par defaut - et on peut definir les parametres des corrections a partir de
fichier ou par argument (normalisation) avec l'appel a des methodes
3/ Gestion de la definition des gains de normalisation globale par argument ou datacard
dans la classe P4AnaParams
4/ Adaptation a la nouvelle interface de P4gvCor ds le constructeur de la classe VisiP4WindowReader
                                                        Reza , 16/05/2018
parent f41e277d
......@@ -13,6 +13,7 @@
#include "p4autils.h"
#include "datacards.h"
#include "strutilxx.h"
#include "xastropack.h"
......@@ -211,7 +212,7 @@ P4AnaParams::P4AnaParams(P4AnaParams const & a)
fgserall_(a.fgserall_), fgtmsel_(a.fgtmsel_), inwsz_(a.inwsz_), fgreorderfreq_(a.fgreorderfreq_),
gain_gnu_file_(a.gain_gnu_file_),
gain_var_file_(a.gain_var_file_), gain_blind_filt_file_(a.gain_blind_filt_file_),
phase_corr_param_file_(a.phase_corr_param_file_),
phase_corr_param_file_(a.phase_corr_param_file_), normalisation_(a.normalisation_),
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_),
......@@ -227,7 +228,7 @@ P4AnaParams& P4AnaParams::operator = (P4AnaParams const & a)
fgserall_=a.fgserall_; fgtmsel_=a. fgtmsel_; inwsz_=a.inwsz_; fgreorderfreq_=a.fgreorderfreq_;
gain_gnu_file_=a.gain_gnu_file_;
gain_var_file_ = a.gain_var_file_; gain_blind_filt_file_ = a.gain_blind_filt_file_;
phase_corr_param_file_ = a.phase_corr_param_file_;
phase_corr_param_file_ = a.phase_corr_param_file_; normalisation_=a.normalisation_;
outfile_=a.outfile_; fitsoutfile_=a.fitsoutfile_;
prtlev_=a.prtlev_; prtmodulo_=a.prtmodulo_;
Nmean_=a.Nmean_; fgTFM_=a.fgTFM_; TFMfreqbin_=a.TFMfreqbin_; TFMtimebin_=a.TFMtimebin_;
......@@ -289,6 +290,16 @@ int P4AnaParams::DecodeArgs(int narg, const char* arg[])
if (narg<2) { cout << "P4AnaParams::DecodeArgs missing/bad argument, -h for help " << endl; return -1; }
phase_corr_param_file_ = arg[2]; arg+=2; narg-=2; lastargs_.clear();
}
else if (fbo=="-norm") {
if (narg<2) { cout << "P4AnaParams::DecodeArgs missing/bad argument, -h for help " << endl; return -1; }
normalisation_.resize((size_t)P4AVisiNumEncoder::getTotalNbVisibilities());
for(size_t jj=0; jj<normalisation_.size(); jj++) normalisation_[jj]=1.;
string sn=arg[2]; vector<string> vs;
SplitStringToVString(sn,vs,',');
size_t jjmx=vs.size(); if (jjmx>normalisation_.size()) jjmx=normalisation_.size();
for(size_t jj=0; jj<jjmx; jj++) normalisation_[jj]=atof(vs[jj].c_str());
arg+=2; narg-=2; lastargs_.clear();
}
else if (fbo=="-o") {
if (narg<2) { cout << "P4AnaParams::DecodeArgs missing/bad argument, -h for help " << endl; return -1; }
outfile_=arg[2]; arg+=2; narg-=2; lastargs_.clear();
......@@ -397,6 +408,13 @@ int P4AnaParams::ReadDatacardFile(const char * dcfilename)
phase_corr_param_file_=dc.SParam("phcor",0,"");
cnt++;
}
if (dc.HasKey("norm")) { // @norm g0 g1 g2 g3 g4 g5 g6 g7
normalisation_.resize((size_t)P4AVisiNumEncoder::getTotalNbVisibilities());
for(size_t jj=0; jj<normalisation_.size(); jj++)
normalisation_[jj]=dc.DParam("norm",jj,1.);
cnt++;
}
if (dc.HasKey("outfile")) { // @outfile OutputFileName
outfile_=dc.SParam("outfile",0,"out.ppf");
cnt++;
......@@ -432,7 +450,8 @@ int P4AnaParams::UsageOptions()
cout<< " P4AnaParams::UsageOptions(): decode datacard and/or command line options \n"
<< " [-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"
<< " [-gnu gain_filename] [-vgf filename] [-bff filename] [-phcor filename] [-norm fac0,fac1...,fac7] \n"
<< " [-o outfilename] [-fo fits_filename] [-nmean N] [-reorderfreq] \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"
......@@ -447,6 +466,7 @@ int P4AnaParams::UsageOptions()
<< " -vgf filename : gain variation parameters file (one line / feed) \n"
<< " -bff filename : blind filetered vectors ppf file (time + value) \n"
<< " -phcor filename : Phase correction parameters filename \n"
<< " -norm fac0,fac1...,fac7 : Define overall normalisation factor for the 8 feeds \n"
<< " -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"
......@@ -483,6 +503,11 @@ ostream& P4AnaParams::Print(ostream& os) const
os<<" Gain, g(nu) file:"<<gain_gnu_file_<<" OutFile="<<outfile_<<endl;
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;
if (normalisation_.size()>0) {
os << " Overall Normalisation[0..7]: ";
for(size_t i=0; i<normalisation_.size(); i++) os << normalisation_[i] << " ; ";
os << 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;
......
......@@ -188,7 +188,9 @@ public:
string gain_blind_filt_file_ ;
//--- nom de fichier de parameteres de corrections des phases
string phase_corr_param_file_;
//--- Facteur de normalisation , 1 valeur / feed , definissable par datacard
vector<double> normalisation_;
//====== Output file definition and print level
//--- nom du fichier de sortie
string outfile_;
......
......@@ -24,44 +24,87 @@
using namespace SOPHYA ;
using namespace std;
P4gvCor::P4gvCor(string const & filename,string const & filt_blind_name)
: parmsfilename_(filename),filt_blind_filename_(filt_blind_name)
P4gvCor::P4gvCor()
: Gt_corr_param_defined_(false), glob_normalisation_defined_(false)
{
InitNormalisation();
}
TArray <double> arr ;
P4gvCor::P4gvCor(string const & param_filename, string const & filt_blind_name)
: Gt_corr_param_defined_(false), glob_normalisation_defined_(false)
{
InitNormalisation();
setGtCorrection(param_filename, filt_blind_name);
}
/*! throws SzMismatchError exception if gnorm.size() != number of feeds (8 for PAON4 */
void P4gvCor::setNormalisationFactors(vector<double> const & gnorm)
{
if (gnorm.size()!=normalisation_.size())
throw SzMismatchError("P4gvCor::setNormalisationFactors(gnorm)/ERROR: Bad gnorm size (!= number of feeds)");
P4AVisiNumEncoder venc;
for(size_t i=0; i<normalisation_.size(); i++) normalisation_[i]=gnorm[i];
for(size_t k=0; k<visi_norm_.size(); k++) {
std::pair<sa_size_t , sa_size_t> pij = venc.getFeedPair(k);
sa_size_t i=pij.first;
sa_size_t j=pij.second;
if (i==j) visi_norm_[k]=normalisation_[i];
else visi_norm_[k]=sqrt(normalisation_[i]*normalisation_[j]);
}
glob_normalisation_defined_=true;
return;
}
ifstream pFile( parmsfilename_.c_str());
if (pFile.is_open())
{
//! Initialises the overall normalisation factor to 1
void P4gvCor::InitNormalisation()
{
P4AVisiNumEncoder venc;
normalisation_.resize(venc.getTotalNbFeeds());
visi_norm_.resize(venc.getTotalNbVisibilities());
vector<double> norms(venc.getTotalNbFeeds());
for(size_t i=0; i<norms.size(); i++) norms[i]=1.;
setNormalisationFactors(norms);
glob_normalisation_defined_=false;
return;
}
sa_size_t nr ;
sa_size_t nc;
arr.ReadASCII (pFile, nr, nc); \
cout << " read "<< nr << " rows with "<< nc << "columns " << endl;
if (nr<8 || nc<3 ) throw IOExc(" P4gvCor : bad number of rows/columns in file ");
for (sa_size_t k=0 ; k < 8 ; k++) {
slopes_.push_back(arr(2,k));
offsets_.push_back(arr(1,k));
cout << " offs "<< arr(1,k) << " sl "<< arr(2,k) << endl;
}
pFile.close();
/*! \param param_filename : name of text file with correction function parameter (one line / feed)
if there are more that 3 parameters on each line (row) , the 4th parameter is considered as the overall normalisation parameter
\param filt_blind_name : filetered blind channel signal PPF filename, signal as a function of time
*/
void P4gvCor::setGtCorrection(string const & param_filename, string const & filt_blind_name)
{
TArray <double> arr ;
ifstream pFile( parmsfilename_.c_str());
if (pFile.is_open()) {
sa_size_t nr ;
sa_size_t nc;
arr.ReadASCII (pFile, nr, nc);
cout << "P4gvCor::setGtCorrection() read "<< nr << " rows with "<< nc << "columns " << endl;
if (nr<8 || nc<3 ) throw IOExc(" P4gvCor::setGtCorrection() : bad number of rows/columns in file ");
vector<double> norm;
for (sa_size_t k=0 ; k < 8 ; k++) {
slopes_.push_back(arr(2,k));
offsets_.push_back(arr(1,k));
cout << k << ": Offset= "<< arr(1,k) << " Slope= "<< arr(2,k);
if (nc>=4) { norm.push_back(arr(3,k)); cout << " Normalisation="<<arr(3,k)<<endl; }
else cout << endl;
}
pFile.close();
}
else {
cout << " P4gvCor : Unable to open file"<<filename << endl;
throw IOExc(" P4gvCor : NO DATA FILE ");
cout << "P4gvCor::setGtCorrection() : Unable to open file"<<param_filename<< endl;
throw IOExc("P4gvCor::setGtCorrection() : NO DATA FILE ");
}
cout << "P4gvCor::setGtCorrection() - reading PPF fileterd blind channel signal from "<<filt_blind_filename_<<endl;
PInPersist pinp(filt_blind_filename_);
pinp >> PPFNameTag("TVFilt_Times") >> filt_vtimes;
pinp >> PPFNameTag("TVFilt_Values") >> filt_values;
vector <double> x;
vector <double> y;
for (sa_size_t k =0 ; k<filt_vtimes.Size() ; k++){
......@@ -69,45 +112,72 @@ P4gvCor::P4gvCor(string const & filename,string const & filt_blind_name)
y.push_back(filt_values(k));
}
my_interp_.DefinePoints (x,y);
}
double P4gvCor::getGain(size_t numch, TimeStamp const & tmstp ){
double P4gvCor::getGain(size_t numch, TimeStamp const & tmstp )
{
if (numch >= slopes_.size() ) throw ParmError(" P4gvCor::GetGain : chanel number outside bounds");
if (!Gt_corr_param_defined_) return 1.;
double vbl = my_interp_(tmstp.ToDays());
double gg = slopes_[numch] * vbl + offsets_[numch];
return(1./gg);
}
void P4gvCor::applyNormalisation(TMatrix <complex <float> > & mtx)
{
if (!glob_normalisation_defined_) return;
if (mtx.NRows() != (sa_size_t)visi_norm_.size()) throw SzMismatchError("P4gvCor::applyNormalisation(mtx) mtx.NRows()!=TotNbVisib");
for(sa_size_t r=0; r< mtx.NRows(); r++) {
mtx.Row(r) *= complex<float>((float)visi_norm_[(size_t)r],0.);
}
return;
}
/*! \param mtx : input/output visibility matrix
\param tmstp : time of the visibility matrix, used to compute the time dependent correction factor G(t)
\param fgdonorm : if true, apply also the overall normalisation factor
*/
void P4gvCor::applyGain(TMatrix <complex <float> > & mtx, TimeStamp const & tmstp, bool fgdonorm)
{
if (!Gt_corr_param_defined_ && (!glob_normalisation_defined_||!fgdonorm)) return;
if (mtx.NRows() != (sa_size_t)visi_norm_.size())
throw SzMismatchError("P4gvCor::applyGain(mtx) mtx.NRows()!=TotNbVisib");
void P4gvCor::applyGain(TMatrix <complex <float> > & mtx, TimeStamp const & tmstp ){
P4AVisiNumEncoder venc;
for(sa_size_t r=0; r< mtx.NRows(); r++) {
std::pair<sa_size_t , sa_size_t> pij = venc.getFeedPair(r);
vector< double > visi_gains(venc.getTotalNbVisibilities());
//--- computing gains
double vbl = my_interp_(tmstp.ToDays());
for(size_t k=0; k<visi_norm_.size(); k++) {
std::pair<sa_size_t , sa_size_t> pij = venc.getFeedPair(k);
sa_size_t i=pij.first;
sa_size_t j=pij.second;
complex <float> gain;
if (i == j ) gain = complex<r_4>( (r_4)getGain (i,tmstp), 0.);
else gain = complex< r_4 >( (r_4) sqrt( getGain(i,tmstp)*getGain(j,tmstp)) , 0.);
cout << " row "<<r<<" = "<< gain <<" "<< getGain (i,tmstp) << " " <<getGain(j,tmstp) << endl;
mtx.Row(r) *= gain ;
double ggi = slopes_[i]*vbl+offsets_[i];
double gg=ggi;
if (i!=j) {
double ggj = slopes_[j]*vbl+offsets_[j];
gg=sqrt(ggi*ggj);
}
if (fgdonorm) gg *= visi_norm_[k];
visi_gains[k]=gg;
}
//--- applying gains ...
for(sa_size_t r=0; r< mtx.NRows(); r++)
mtx.Row(r) *= complex<float>((float)visi_gains[(size_t)r],0.);
return;
}
}
......@@ -24,13 +24,33 @@ using namespace std;
//---------------------------------------------------------------------------------
//---- Gestion de la cor des variations de gain avec la voie 4V pour PAON-4
//! \brief class for managing relative gain time variations with 4V channel for PAON-4
/*!
\brief class for managing time variations of the overall gain G(t) for PAON-4
A blind channel (4V for some period) is used to track temperature dependent gain
*/
class P4gvCor {
public:
P4gvCor(string const & filename,string const & filt_blind_name);
//! Default constructor
P4gvCor();
/*! \brief Constructor, with the specification of two files for the corrections with time/temperature
\sa setGtCorrection */
P4gvCor(string const & param_filename,string const & filt_blind_name);
//! Define the overall normalisation factors for the auto-correlation channels
void setNormalisationFactors(vector<double> const & gnorm);
//! Define the time dependent gain correction (correction function parameters, and filetered blind channel filename
void setGtCorrection(string const & param_filename,string const & filt_blind_name);
//! Return the time/temperature dependent gain factor for channel numch and time tmstp
double getGain(size_t numch, TimeStamp const & tmstp );
void applyGain(TMatrix <complex <float> > & mtx, TimeStamp const & tmstp );
//! apply the overall normalisation factors to the visibility matrix
void applyNormalisation(TMatrix <complex <float> > & mtx);
//! apply the time/temperature dependent gain factor and normalisation factor to the visibility matrix at time tmstp and optionaly normalisation
void applyGain(TMatrix <complex <float> > & mtx, TimeStamp const & tmstp, bool fgdonorm=true);
protected:
void InitNormalisation();
//--- variables pour la correction du gain dependant du temps / temperature G(t)
bool Gt_corr_param_defined_; // true : parametres pour la correction de gain G(t) defini
std::string parmsfilename_;
std::string filt_blind_filename_;
vector <double> slopes_ ;
......@@ -39,7 +59,10 @@ protected:
TVector<double> filt_vtimes;
TVector<double> filt_values;
SLinInterp1D my_interp_;
//--- Variables pour l'application de la normalisation globale
bool glob_normalisation_defined_; // true : facteur de normalisation global a ete defini
vector <double> normalisation_; // facteur de normalisation (1fac / voie=autocorrelation)
vector < double > visi_norm_; // facteur de normalisation (1fac / visibilite = paire de voies)
};
......
......@@ -56,8 +56,11 @@ VisiP4WindowReader::VisiP4WindowReader(P4AnaParams & params)
p4gnu_= new P4gnuGain(params.gain_gnu_file_);
fgown_p4gnu_=true;
}
if ( params.gain_var_file_.length()>1 && params.gain_blind_filt_file_.length() > 1 ) {
p4gvf_ = new P4gvCor(params.gain_var_file_, params.gain_blind_filt_file_);
if ( ((params.gain_var_file_.length()>1)&&(params.gain_blind_filt_file_.length()>1)) || (params.normalisation_.size()>0) ) {
p4gvf_ = new P4gvCor;
if ((params.gain_var_file_.length()>1)&&(params.gain_blind_filt_file_.length()>1))
p4gvf_->setGtCorrection(params.gain_var_file_, params.gain_blind_filt_file_);
if (params.normalisation_.size()>0) p4gvf_->setNormalisationFactors(params.normalisation_);
fgown_p4gvf_=true;
}
if (params.phase_corr_param_file_.length() > 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