Commit 4a0fbe47 authored by Reza  ANSARI's avatar Reza ANSARI
Browse files

Ajout de l'option -dosigma, pour rendre le calcul des cartes TFM...

Ajout de l'option -dosigma, pour rendre le calcul des cartes TFM (Time-Frequency) de la dispersion (sigma) des differents signaux et sa prise en charge ds visiavg.cc , Reza 15/07/2018
parent 5452b8e7
......@@ -40,3 +40,5 @@ norm g1H g2H g3H g4H g1V g2V g3V g4V
### Select VV and/or HV cross correlation processing
@dovv
@dohv
### Require computation of dispersion (sigma) TFM maps, in addition to mean (average) TFM maps
@dosigma
......@@ -214,7 +214,8 @@ P4AnaParams::P4AnaParams()
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),
fgRA_(false), RAcenter_(12.), RAwidth_(24.), RAresol_mn_(1.), doVV_(false),doHV_(false)
fgRA_(false), RAcenter_(12.), RAwidth_(24.), RAresol_mn_(1.),
doVV_(false), doHV_(false), doSigma_(false)
{
}
......@@ -228,7 +229,8 @@ P4AnaParams::P4AnaParams(P4AnaParams const & a)
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_),
fgRA_(a.fgRA_), RAcenter_(a.RAcenter_), RAwidth_(a.RAwidth_), RAresol_mn_(a.RAresol_mn_),doVV_(a.doVV_),doHV_(a.doHV_),
fgRA_(a.fgRA_), RAcenter_(a.RAcenter_), RAwidth_(a.RAwidth_), RAresol_mn_(a.RAresol_mn_),
doVV_(a.doVV_), doHV_(a.doHV_), doSigma_(a.doSigma_),
lastargs_(a.lastargs_)
{
}
......@@ -245,7 +247,8 @@ P4AnaParams& P4AnaParams::operator = (P4AnaParams const & a)
prtlev_=a.prtlev_; prtmodulo_=a.prtmodulo_;
Nmean_=a.Nmean_; fgTFM_=a.fgTFM_; TFMfreqbin_=a.TFMfreqbin_; TFMtimebin_=a.TFMtimebin_;
fgdtable_=a.fgdtable_; fbands_=a.fbands_;
fgRA_=a.fgRA_; RAcenter_=a.RAcenter_; RAwidth_=a.RAwidth_; RAresol_mn_=a.RAresol_mn_;doVV_=a.doVV_;doHV_=a.doHV_;
fgRA_=a.fgRA_; RAcenter_=a.RAcenter_; RAwidth_=a.RAwidth_; RAresol_mn_=a.RAresol_mn_;
doVV_=a.doVV_; doHV_=a.doHV_; doSigma_=a.doSigma_;
lastargs_=a.lastargs_;
return *this;
}
......@@ -290,7 +293,12 @@ int P4AnaParams::DecodeArgs(int narg, const char* arg[])
doVV_=true; arg++; narg--; lastargs_.clear();
}
else if (fbo=="-dohv") {
cout << "P4AnaParams::DecodeArgs : HV cross correlations demanded - may need lof of memory ! "<<endl;
doHV_=true; arg++; narg--; lastargs_.clear();
}
else if (fbo=="-dosigma") {
doSigma_=true; arg++; narg--; lastargs_.clear();
}
else if (fbo=="-dohv") {
doHV_=true; arg++; narg--; lastargs_.clear();
}
else if (fbo=="-gnu") {
......@@ -458,13 +466,15 @@ int P4AnaParams::ReadDatacardFile(const char * dcfilename)
prtlev_=dc.IParam("prtlev",0,1); prtmodulo_=dc.IParam("prtlev",1,1);
cnt++;
}
if (dc.HasKey("dovv")) { //@dovv
if (dc.HasKey("dovv")) { //@dovv : process also VV cross-correlations
doVV_=true; cnt++;
}
if (dc.HasKey("dohv")) { //@dohv
if (dc.HasKey("dohv")) { //@dohv : process also HV cross-correlations
doHV_=true; cnt++;
}
if (dc.HasKey("dosigma")) { //@dosigma : computes also dispersion (Sigma) TFM maps
doSigma_=true; cnt++;
}
for (size_t i=0; i<dc_fbands_.size(); i++) fbands_.push_back(dc_fbands_[i]);
cout << " P4AnaParams::ReadDatacardFile(FileName=" << dcfilename << " ) -> decoded " << cnt+dc_fbands_.size()<<" cards"<<endl;
dc_fbands_.clear();
......@@ -479,7 +489,7 @@ int P4AnaParams::UsageOptions()
<< " [-gnu gain_filename] [-vgf filename] [-bff filename] [-phcor filename] [-norm fac1H,fac2H...,fac4V] \n"
<< " [-o outfilename] [-fo fits_filename] [-nmean N] [-reorderfreq] \n"
<< " [-tfm timebin,freqbin] [-ra RAcenterHours,RAwidthHours,RAresolMinutes] [-freq f0,df] \n"
<< " [-dovv] [-dohv] [-prt level] \n"
<< " [-dovv] [-dohv] [-dosigma] [-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"
......@@ -503,6 +513,7 @@ int P4AnaParams::UsageOptions()
<< " several frequency bands can be specified by mutiple -freq options \n"
<< " -dovv : processing also VV cross correlations \n"
<< " -dohv : processing also HV cross correlations \n"
<< " -dosigma : compute also TFM maps of dispersion (sigma) \n"
<< " -prt level,modulo: define print level and print modulo count \n"
<< endl;
return 0;
......@@ -545,7 +556,7 @@ ostream& P4AnaParams::Print(ostream& os) const
else for(size_t i=0; i<fbands_.size(); i++) os<<"FreqBand["<<i<<"]: "<<fbands_[i]<<endl;
if (doVV_) cout << "* Processing VV cross correlations demanded "<<endl;
if (doHV_) cout << "* Processing HV cross correlations demanded "<<endl;
if (doSigma_) cout << "* Computing also dispersion (Sigma) TFM maps, in addition to mean (average) "<<endl;
return os;
}
......
......@@ -218,6 +218,7 @@ public:
bool doVV_; // true : also computes VxV cross correlations tf maps (visiavg)
bool doHV_; // true : also computes HxV cross correlations tf maps (visiavg)
bool doSigma_; // true : computes also dispersion (sigma) tf maps , in addition to mean (visiavg)
};
inline ostream& operator << (ostream& os, P4AnaParams const& a)
......
......@@ -195,24 +195,30 @@ int main(int narg, const char* arg[])
if (cnt==0) { //resizing matrices for sum of auto-correlations and sum of 6 cross-correlations
acsum.SetSize(8, vismtx.NCols());
acsum_sq.SetSize(8, vismtx.NCols());
if (params.doSigma_) // for computing Sigma, if required
acsum_sq.SetSize(8, vismtx.NCols());
cxsum.SetSize(6, vismtx.NCols());
cxsum_sq_rp.SetSize(6, vismtx.NCols());
cxsum_sq_ip.SetSize(6, vismtx.NCols());
if (params.doSigma_) { // for computing Sigma, if required
cxsum_sq_rp.SetSize(6, vismtx.NCols());
cxsum_sq_ip.SetSize(6, vismtx.NCols());
}
// VV if needed
if (params.doVV_){
cxsum_vv.SetSize(6, vismtx.NCols());
cxsum_vv_sq_rp.SetSize(6, vismtx.NCols());
cxsum_vv_sq_ip.SetSize(6, vismtx.NCols());
if (params.doSigma_) { // for computing Sigma, if required
cxsum_vv_sq_rp.SetSize(6, vismtx.NCols());
cxsum_vv_sq_ip.SetSize(6, vismtx.NCols());
}
}
// HV if needed
if (params.doHV_){
cxsum_hv.SetSize(16, vismtx.NCols());
cxsum_hv_sq_rp.SetSize(16, vismtx.NCols());
cxsum_hv_sq_ip.SetSize(16, vismtx.NCols());
if (params.doSigma_) { // for computing Sigma, if required
cxsum_hv_sq_rp.SetSize(16, vismtx.NCols());
cxsum_hv_sq_ip.SetSize(16, vismtx.NCols());
}
}
......@@ -220,24 +226,29 @@ int main(int narg, const char* arg[])
tfmSY=vismtx.NCols()/TFMfbin;
//allocating 8 Auto-Corr time-frequency maps
cout<<"visiavg/Info: allocating 8 AutoCor Time-Frequency maps : Time->NX="<<tfmSX<<" x Freq->NY="<<tfmSY<<endl;
for(int k=0; k<8; k++) vtfmac.push_back( TArray< r_4 >(tfmSX, tfmSY) );
cout <<" and 8 for the the variance maps "<<endl;
for(int k=0; k<8; k++) vtfmac_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
for(int k=0; k<8; k++) vtfmac.push_back( TArray< r_4 >(tfmSX, tfmSY) );
if (params.doSigma_) { // for computing Sigma, if required
cout <<" and 8 for the the variance maps "<<endl;
for(int k=0; k<8; k++) vtfmac_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
}
//allocating 6 Cross-Corr H-H time-frequency maps
cout<<"visiavg/Info: allocating H-H cross-cor Time-Frequency maps : Time->NX="<<tfmSX<<" x Freq->NY="<<tfmSY<<endl;
for(int k=0; k<6; k++) vtfm.push_back( TArray< complex<r_4> >(tfmSX, tfmSY) );
cout << "and the 2x6 for squares of real & imaginary parts " << endl;
for(int k=0; k<6; k++) vtfm_rp_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
for(int k=0; k<6; k++) vtfm_ip_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
if (params.doSigma_) { // for computing Sigma, if required
cout << "and the 2x6 for squares of real & imaginary parts " << endl;
for(int k=0; k<6; k++) vtfm_rp_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
for(int k=0; k<6; k++) vtfm_ip_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
}
// VV if needed
if (params.doVV_){
//allocating 6 Cross-Corr V-V time-frequency maps
cout<<"visiavg/Info: allocating V-V cross-cor Time-Frequency maps : Time->NX="<<tfmSX<<" x Freq->NY="<<tfmSY<<endl;
for(int k=0; k<6; k++) vtfm_vv.push_back( TArray< complex<r_4> >(tfmSX, tfmSY) );
cout << "and the 2x6 for squares of real & imaginary parts " << endl;
for(int k=0; k<6; k++) vtfm_vv_rp_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
for(int k=0; k<6; k++) vtfm_vv_ip_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
if (params.doSigma_) { // for computing Sigma, if required
cout << "and the 2x6 for squares of real & imaginary parts " << endl;
for(int k=0; k<6; k++) vtfm_vv_rp_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
for(int k=0; k<6; k++) vtfm_vv_ip_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
}
}
// HV if needed
......@@ -245,9 +256,11 @@ int main(int narg, const char* arg[])
//allocating 16 Cross-Corr H-V time-frequency maps
cout<<"visiavg/Info: allocating 16 H-V cross-cor Time-Frequency maps : Time->NX="<<tfmSX<<" x Freq->NY="<<tfmSY<<endl;
for(int k=0; k<16; k++) vtfm_hv.push_back( TArray< complex<r_4> >(tfmSX, tfmSY) );
cout << "and the 2x16 for squares of real & imaginary parts " << endl;
for(int k=0; k<16; k++) vtfm_hv_rp_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
for(int k=0; k<16; k++) vtfm_hv_ip_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
if (params.doSigma_) { // for computing Sigma, if required
cout << "and the 2x16 for squares of real & imaginary parts " << endl;
for(int k=0; k<16; k++) vtfm_hv_rp_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
for(int k=0; k<16; k++) vtfm_hv_ip_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
}
}
......@@ -260,147 +273,159 @@ int main(int narg, const char* arg[])
if (prtlev>0)
cout<<"visiavg/Info: dateobs="<<dateobs<<" SecondsPart()="<<dateobs.SecondsPart()<<endl;
acsum = complex<r_4>(0.,0.);
acsum_sq = 0.;
cxsum = complex<r_4>(0.,0.);
cxsum_sq_rp = 0.;
cxsum_sq_ip = 0.;
if (params.doSigma_) { // for computing Sigma, if required
acsum_sq = 0.;
cxsum_sq_rp = 0.;
cxsum_sq_ip = 0.;
}
// VV if needed
if (params.doVV_){
cxsum_vv = complex<r_4>(0.,0.);
cxsum_vv_sq_rp = 0.;
cxsum_vv_sq_ip = 0.;
if (params.doSigma_) { // for computing Sigma, if required
cxsum_vv_sq_rp = 0.;
cxsum_vv_sq_ip = 0.;
}
}
// HV if needed
if (params.doHV_){
cxsum_hv = complex<r_4>(0.,0.);
cxsum_hv_sq_rp = 0.;
cxsum_hv_sq_ip = 0.;
if (params.doSigma_) { // for computing Sigma, if required
cxsum_hv_sq_rp = 0.;
cxsum_hv_sq_ip = 0.;
}
}
}
// sum (integration) along the time axis
for(size_t k=0; k<KVAC.size(); k++) acsum.Row(k) += vismtx.Row(KVAC[k]); // Les auto-correlations
for(size_t k=0; k<KVAC.size(); k++) {
TVector<r_4> tmp = real(vismtx.Row(KVAC[k]));
acsum_sq.Row(k) += tmp.MulElt(tmp,tmp) ;
if (params.doSigma_) { // for computing Sigma, if required
for(size_t k=0; k<KVAC.size(); k++) {
TVector<r_4> tmp = real(vismtx.Row(KVAC[k]));
acsum_sq.Row(k) += tmp.MulElt(tmp,tmp) ;
}
}
for(size_t k=0; k<KVCXHH.size(); k++){
cxsum.Row(k) += vismtx.Row(KVCXHH[k]); // les cross-correlations
TVector<r_4> tmp = real(vismtx.Row(KVCXHH[k]));
cxsum_sq_rp.Row(k) += tmp.MulElt(tmp,tmp) ;
tmp = imag(vismtx.Row(KVCXHH[k]));
cxsum_sq_ip.Row(k) += tmp.MulElt(tmp,tmp) ;
if (params.doSigma_) { // for computing Sigma, if required
TVector<r_4> tmp = real(vismtx.Row(KVCXHH[k]));
cxsum_sq_rp.Row(k) += tmp.MulElt(tmp,tmp) ;
tmp = imag(vismtx.Row(KVCXHH[k]));
cxsum_sq_ip.Row(k) += tmp.MulElt(tmp,tmp) ;
}
}
if (params.doVV_){
for(size_t k=0; k<KVCXVV.size(); k++){
cxsum_vv.Row(k) += vismtx.Row(KVCXVV[k]); // les cross-correlations
TVector<r_4> tmp = real(vismtx.Row(KVCXVV[k]));
cxsum_vv_sq_rp.Row(k) += tmp.MulElt(tmp,tmp) ;
tmp = imag(vismtx.Row(KVCXVV[k]));
cxsum_vv_sq_ip.Row(k) += tmp.MulElt(tmp,tmp) ;
if (params.doSigma_) { // for computing Sigma, if required
TVector<r_4> tmp = real(vismtx.Row(KVCXVV[k]));
cxsum_vv_sq_rp.Row(k) += tmp.MulElt(tmp,tmp) ;
tmp = imag(vismtx.Row(KVCXVV[k]));
cxsum_vv_sq_ip.Row(k) += tmp.MulElt(tmp,tmp) ;
}
}
}
if (params.doHV_){
for(size_t k=0; k<KVCXHV.size(); k++){
cxsum_hv.Row(k) += vismtx.Row(KVCXHV[k]); // les cross-correlations HV
TVector<r_4> tmp = real(vismtx.Row(KVCXHV[k]));
cxsum_hv_sq_rp.Row(k) += tmp.MulElt(tmp,tmp) ;
tmp = imag(vismtx.Row(KVCXHV[k]));
cxsum_hv_sq_ip.Row(k) += tmp.MulElt(tmp,tmp) ;
if (params.doSigma_) { // for computing Sigma, if required
TVector<r_4> tmp = real(vismtx.Row(KVCXHV[k]));
cxsum_hv_sq_rp.Row(k) += tmp.MulElt(tmp,tmp) ;
tmp = imag(vismtx.Row(KVCXHV[k]));
cxsum_hv_sq_ip.Row(k) += tmp.MulElt(tmp,tmp) ;
}
}
}
I++; // incrementing DeltaTime counter
if (I==deltaIavg) {
// 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 << "visiavg/Warning: something wrong in the logic , (TFMtmidx="<<TFMtmidx<<") >= (tfmSX="<<tfmSX<<")"
<< " for read count="<<cnt<<endl;
}
else if (I==deltaIavg) { // Filling TFM maps
//---- 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(acsum.Row(k));
TVector<r_4> vacsq = acsum_sq.Row(k);
if (TFMtmidx<tfmSX) { // we check that our time index did not go beyond the allocated array size (might not be necessary)
TArray< r_4 > & tfmap = vtfmac[k];
TArray< r_4 > & tfmap = vtfmac[k];
for(sa_size_t jy=0; jy<tfmSY; jy++) { // frequency binning
tfmap(TFMtmidx, jy) = vac( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
}
if (params.doSigma_) { // for computing Sigma, if required
TVector<r_4> vacsq = acsum_sq.Row(k);
TArray< r_4 > & tfmap_sq = vtfmac_sq[k];
for(sa_size_t jy=0; jy<tfmSY; jy++) { // frequency binning
tfmap(TFMtmidx, jy) = vac( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
tfmap_sq(TFMtmidx, jy) = vacsq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
}
}
}
} //-- end of computing sigmas
} //----- 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 = cxsum.Row(k);
TVector<r_4> vcxprsq = cxsum_sq_rp.Row(k);
TVector<r_4> vcxpisq = cxsum_sq_ip.Row(k);
if (TFMtmidx<tfmSX) { // we check that our time index did not go beyond the allocated array size (might not be necessary)
TArray< complex<r_4> > & tfmap = vtfm[k];
TArray< complex<r_4> > & tfmap = vtfm[k];
for(sa_size_t jy=0; jy<tfmSY; jy++) {
tfmap(TFMtmidx, jy) = vcx( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
}
if (params.doSigma_) { // for computing Sigma, if required
TVector<r_4> vcxprsq = cxsum_sq_rp.Row(k);
TVector<r_4> vcxpisq = cxsum_sq_ip.Row(k);
TArray< r_4 > & tfmapsqpr = vtfm_rp_sq[k];
TArray< r_4 > & tfmapsqpi = vtfm_ip_sq[k];
for(sa_size_t jy=0; jy<tfmSY; jy++) {
tfmap(TFMtmidx, jy) = vcx( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
tfmapsqpr(TFMtmidx, jy) = vcxprsq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
tfmapsqpi(TFMtmidx, jy) = vcxpisq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
}
}
} //-- end of computing sigmas
} //----- end of loop over the 6 Xcor
if (params.doVV_){ // Option VV
//---- On s'occupe des 6 cross-correlations 1V-2V ... 3V-4V
for(size_t k=0; k<KVCXVV.size(); k++) { // loop over the 6 Xcor
TVector< complex<r_4> > vcx = cxsum_vv.Row(k);
TVector<r_4> vcxprsq = cxsum_vv_sq_rp.Row(k);
TVector<r_4> vcxpisq = cxsum_vv_sq_ip.Row(k);
if (TFMtmidx<tfmSX) { // we check that our time index did not go beyond the allocated array size (might not be necessary)
TArray< complex<r_4> > & tfmap = vtfm_vv[k];
TArray< complex<r_4> > & tfmap = vtfm_vv[k];
for(sa_size_t jy=0; jy<tfmSY; jy++) {
tfmap(TFMtmidx, jy) = vcx( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
}
if (params.doSigma_) { // for computing Sigma, if required
TVector<r_4> vcxprsq = cxsum_vv_sq_rp.Row(k);
TVector<r_4> vcxpisq = cxsum_vv_sq_ip.Row(k);
TArray< r_4 > & tfmapsqpr = vtfm_vv_rp_sq[k];
TArray< r_4 > & tfmapsqpi = vtfm_vv_ip_sq[k];
for(sa_size_t jy=0; jy<tfmSY; jy++) {
tfmap(TFMtmidx, jy) = vcx( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
tfmapsqpr(TFMtmidx, jy) = vcxprsq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
tfmapsqpi(TFMtmidx, jy) = vcxpisq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
}
}
}
} //-- end of computing sigmas
} //----- end of loop over the 6 Xcor VV
} // end VV option
if (params.doHV_){ // Option HV
//---- On s'occupe des 16 cross-correlations 1H-1V ... 4H-4V
for(size_t k=0; k<KVCXHV.size(); k++) { // loop over the 16 Xcor
TVector< complex<r_4> > vcx = cxsum_hv.Row(k);
TVector<r_4> vcxprsq = cxsum_hv_sq_rp.Row(k);
TVector<r_4> vcxpisq = cxsum_hv_sq_ip.Row(k);
if (TFMtmidx<tfmSX) { // we check that our time index did not go beyond the allocated array size (might not be necessary)
TArray< complex<r_4> > & tfmap = vtfm_hv[k];
TArray< complex<r_4> > & tfmap = vtfm_hv[k];
for(sa_size_t jy=0; jy<tfmSY; jy++) {
tfmap(TFMtmidx, jy) = vcx( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
}
if (params.doSigma_) { // for computing Sigma, if required
TVector<r_4> vcxprsq = cxsum_hv_sq_rp.Row(k);
TVector<r_4> vcxpisq = cxsum_hv_sq_ip.Row(k);
TArray< r_4 > & tfmapsqpr = vtfm_hv_rp_sq[k];
TArray< r_4 > & tfmapsqpi = vtfm_hv_ip_sq[k];
for(sa_size_t jy=0; jy<tfmSY; jy++) {
tfmap(TFMtmidx, jy) = vcx( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
tfmapsqpr(TFMtmidx, jy) = vcxprsq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
tfmapsqpi(TFMtmidx, jy) = vcxpisq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
}
}
} //-- end of computing sigmas
} //----- end of loop over the 16 Xcor HV
} // end HV option
double tdif = cfdate.TimeDifferenceSeconds(cfdate,dateobs)/2.;
......@@ -412,11 +437,12 @@ int main(int narg, const char* arg[])
}
cnt++;
if ((cnt>0)&&(cntnt%10==0)&&(cntnt>pcntnt)) {
cout<<"visiavg/Info: TFM-Map fill cnt="<<cntnt<<" VisMtxCount="<<cnt<<" /Max="<<Imax<<" DateObs="<<dateobs<<endl;
cout<<"visiavg/Info: TFM-Map fill cnt="<<cntnt<<" VisMtxCount="<<cnt
<<" /Max="<<wreader.getTotalNbWindows()<<" DateObs="<<dateobs<<endl;
pcntnt=cntnt;
}
}
cout<<"visiavg/Info: count="<<cnt<<" visimtx read "<<endl;
cout<<"visiavg/Info: count="<<cnt*wreader.getWindowSize()<<" Visibility Matrices read "<<endl;
// --- Sauvegarde cartes temps-frequence en fits
//FitsABTWriter * fbtw = NULL;
......@@ -444,15 +470,16 @@ int main(int narg, const char* arg[])
ext_names.push_back(tfm_names[k]);
(*fos)<< tfmap;
}
TArray< r_4 > & tfmapsq = vtfmac_sq[k];
tfmapsq *= (r_4)(1./((double)deltaIavg*(double)TFMfbin));
tfmapsq = tfmapsq - tfmap.MulElt(tfmap,tfmap) ; // tfmap ->tfmap*tfmap !!
potfm << PPFNameTag(tfmsq_names[k]) << tfmapsq;
if (fos != NULL) {
ext_names.push_back(tfmsq_names[k]);
(*fos)<< tfmapsq;
}
if (params.doSigma_) { // for saving TFM-Sigma, if required
TArray< r_4 > & tfmapsq = vtfmac_sq[k];
tfmapsq *= (r_4)(1./((double)deltaIavg*(double)TFMfbin));
tfmapsq = tfmapsq - tfmap.MulElt(tfmap,tfmap) ; // tfmap ->tfmap*tfmap !!
potfm << PPFNameTag(tfmsq_names[k]) << tfmapsq;
if (fos != NULL) {
ext_names.push_back(tfmsq_names[k]);
(*fos)<< tfmapsq;
}
} // end of saving TFM-of-sigma
}
// --- renormalizing and saving H-H Cross-Corr time-frequency maps
......@@ -475,24 +502,25 @@ int main(int narg, const char* arg[])
ext_names.push_back(string(tfmCC_names[k])+"_imag");
(*fos)<< imag(tfmap);
}
TArray< r_4 > tfmapr = real(tfmap);
TArray< r_4 > tfmapi = imag(tfmap);
tfmap_sqpr *= ((r_4)(1./((double)deltaIavg*(double)TFMfbin)));
tfmap_sqpi *= ((r_4)(1./((double)deltaIavg*(double)TFMfbin)));
tfmapr = tfmapr.MulElt(tfmapr,tfmapr) ;
tfmap_sqpr -= tfmapr ;
tfmapi = tfmapi.MulElt(tfmapi,tfmapi) ;
tfmap_sqpi -= tfmapi;
potfm << PPFNameTag(vrtfmCC_names[k]) << tfmap_sqpr;
potfm << PPFNameTag(vitfmCC_names[k]) << tfmap_sqpi;
if (fos != NULL) {
ext_names.push_back(vrtfmCC_names[k]);
(*fos)<< tfmap_sqpr;
ext_names.push_back(vitfmCC_names[k]);
(*fos)<< tfmap_sqpi;
}
if (params.doSigma_) { // for saving TFM-Sigma, if required
TArray< r_4 > tfmapr = real(tfmap);
TArray< r_4 > tfmapi = imag(tfmap);
tfmap_sqpr *= ((r_4)(1./((double)deltaIavg*(double)TFMfbin)));
tfmap_sqpi *= ((r_4)(1./((double)deltaIavg*(double)TFMfbin)));
tfmapr = tfmapr.MulElt(tfmapr,tfmapr) ;
tfmap_sqpr -= tfmapr ;
tfmapi = tfmapi.MulElt(tfmapi,tfmapi) ;
tfmap_sqpi -= tfmapi;
potfm << PPFNameTag(vrtfmCC_names[k]) << tfmap_sqpr;
potfm << PPFNameTag(vitfmCC_names[k]) << tfmap_sqpi;
if (fos != NULL) {
ext_names.push_back(vrtfmCC_names[k]);
(*fos)<< tfmap_sqpr;
ext_names.push_back(vitfmCC_names[k]);
(*fos)<< tfmap_sqpi;
}
} // end of saving TFM-of-sigma
}
// ========================== end saving HH XCorr
......@@ -519,26 +547,26 @@ int main(int narg, const char* arg[])
ext_names.push_back(string(tfmVV_names[k])+"_imag");
(*fos)<< imag(tfmap);
}
TArray< r_4 > tfmapr = real(tfmap);
TArray< r_4 > tfmapi = imag(tfmap);
tfmap_sqpr *= ((r_4)(1./((double)deltaIavg*(double)TFMfbin)));
tfmap_sqpi *= ((r_4)(1./((double)deltaIavg*(double)TFMfbin)));
tfmapr = tfmapr.MulElt(tfmapr,tfmapr) ;
tfmap_sqpr -= tfmapr ;
tfmapi = tfmapi.MulElt(tfmapi,tfmapi) ;
tfmap_sqpi -= tfmapi;
potfm << PPFNameTag(vrtfmVV_names[k]) << tfmap_sqpr;
potfm << PPFNameTag(vitfmVV_names[k]) << tfmap_sqpi;
if (fos != NULL) {
ext_names.push_back(vrtfmVV_names[k]);
(*fos)<< tfmap_sqpr;
ext_names.push_back(vitfmVV_names[k]);
(*fos)<< tfmap_sqpi;
}
if (params.doSigma_) { // for saving TFM-Sigma, if required
TArray< r_4 > tfmapr = real(tfmap);
TArray< r_4 > tfmapi = imag(tfmap);
tfmap_sqpr *= ((r_4)(1./((double)deltaIavg*(double)TFMfbin)));
tfmap_sqpi *= ((r_4)(1./((double)deltaIavg*(double)TFMfbin)));
tfmapr = tfmapr.MulElt(tfmapr,tfmapr) ;
tfmap_sqpr -= tfmapr ;
tfmapi = tfmapi.MulElt(tfmapi,tfmapi) ;
tfmap_sqpi -= tfmapi;
potfm << PPFNameTag(vrtfmVV_names[k]) << tfmap_sqpr;
potfm << PPFNameTag(vitfmVV_names[k]) << tfmap_sqpi;
if (fos != NULL) {
ext_names.push_back(vrtfmVV_names[k]);
(*fos)<< tfmap_sqpr;
ext_names.push_back(vitfmVV_names[k]);
(*fos)<< tfmap_sqpi;
}
} // end of saving TFM-of-sigma
}
} // end option VV
......@@ -577,24 +605,26 @@ int main(int narg, const char* arg[])
ext_names.push_back(string(tfmHV_names[k])+"_imag");
(*fos)<< imag(tfmap);
}
TArray< r_4 > tfmapr = real(tfmap);
TArray< r_4 > tfmapi = imag(tfmap);
tfmap_sqpr *= ((r_4)(1./((double)deltaIavg*(double)TFMfbin)));
tfmap_sqpi *= ((r_4)(1./((double)deltaIavg*(double)TFMfbin)));
tfmapr = tfmapr.MulElt(tfmapr,tfmapr) ;
tfmap_sqpr -= tfmapr ;
tfmapi = tfmapi.MulElt(tfmapi,tfmapi) ;
tfmap_sqpi -= tfmapi;
potfm << PPFNameTag(vrtfmHV_names[k]) << tfmap_sqpr;
potfm << PPFNameTag(vitfmHV_names[k]) << tfmap_sqpi;
if (fos != NULL) {
ext_names.push_back(vrtfmHV_names[k]);
(*fos)<< tfmap_sqpr;
ext_names.push_back(vitfmHV_names[k]);
(*fos)<< tfmap_sqpi;
}
if (params.doSigma_) { // for saving TFM-Sigma, if required
TArray< r_4 > tfmapr = real(tfmap);
TArray< r_4 > tfmapi = imag(tfmap);
tfmap_sqpr *= ((r_4)(1./((double)deltaIavg*(double)TFMfbin)));
tfmap_sqpi *= ((r_4)(1./((double)deltaIavg*(double)TFMfbin)));
tfmapr = tfmapr.MulElt(tfmapr,tfmapr) ;
tfmap_sqpr -= tfmapr ;
tfmapi = tfmapi.MulElt(tfmapi,tfmapi) ;
tfmap_sqpi -= tfmapi;
potfm << PPFNameTag(vrtfmHV_names[k]) << tfmap_sqpr;
potfm << PPFNameTag(vitfmHV_names[k]) << tfmap_sqpi;
if (fos != NULL) {
ext_names.push_back(vrtfmHV_names[k]);
(*fos)<< tfmap_sqpr;
ext_names.push_back(vitfmHV_names[k]);
(*fos)<< tfmap_sqpi;
}
} // end of saving TFM-of-sigma
}
} // end option HV
......