visiavg.cc 16.2 KB
Newer Older
OP's avatar
OP committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
// Utilisation de SOPHYA pour faciliter les tests ...
#include "sopnamsp.h"
#include "machdefs.h"

/* ---------------------------------------------------------- 
   Projet BAORadio/PAON4 - (C) LAL/IRFU  2017

   visiavg: programme de lecture des fichiers matrices de 
   visibilites de PAON4, calcul de visibilities moyennes 
    en bin de temps et de frequence
   O. Perdereau, R.Ansari   -  LAL
   ---------------------------------------------------------- */

// 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 "p4autils.h"
41
#include "visp4winreader.h"
OP's avatar
OP committed
42
#include "p4gnugain.h"
OP's avatar
modifs  
OP committed
43
#include "p4gvcor.h"
44
#include "p4phcor.h"
OP's avatar
OP committed
45

OP's avatar
OP committed
46 47 48
#include "fitsioserver.h"
#include "fiosinit.h"

OP's avatar
OP committed
49 50 51 52 53 54 55 56 57 58 59 60 61 62 63
int Usage(void);
int Usage(void)
{
  cout << " --- visiavg.cc : Read PPF files produced by mfacq time-frequency\n" << endl;
  cout << " Usage: visiavg [-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();
OP's avatar
OP committed
64
  FitsIOServerInit();
OP's avatar
OP committed
65 66
  P4AnaParams params;
  params.DecodeArgs(narg, arg);
OP's avatar
OP committed
67

OP's avatar
OP committed
68 69
  string outfile = params.outfile_;
  if (outfile.length()<1)  outfile = "visavg.ppf";
OP's avatar
OP committed
70

OP's avatar
OP committed
71 72 73 74
  string fitsoutfile = params.fitsoutfile_;
  if (fitsoutfile.length()>=1) {
    fitsoutfile = "!"+fitsoutfile ; // adds '!' ?
  }
OP's avatar
modifs  
OP committed
75

OP's avatar
OP committed
76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106
  int deltaIavg = params.TFMtimebin_;
  sa_size_t TFMfbin = params.TFMfreqbin_;
  int Imin = params.Imin_, Imax = params.Imax_, Istep = params.Istep_; 
  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 <<"visiavg/Info: Path BAO5:"<<params.inpath5_<<" BAO6:"<<params.inpath6_<<"\n"
       <<"fgreorderfreq="<<params.fgreorderfreq_<<"\n"
       <<"Imin,max,step="<<Imin<<","<<Imax<<","<<Istep<<" DeltaIAvg="<<deltaIavg<<"\n"
       <<"outfile="<<outfile<<" PrtLev="<<prtlev<<endl;

  if (!FgTFM) {
    cout<<" visiavg/parameter error : specify Time-Frequency map parameter with -tfm "<<endl;
    return 5;
  }

  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;
  }
OP's avatar
modifs  
OP committed
107 108 109 110 111 112 113 114 115

  string gvfile = params.gain_var_file_;
  string fltfile = params.gain_blind_filt_file_;
  P4gvCor * p4gvf = NULL ;
  cout <<" gvfile "<< gvfile << endl;
  cout <<"fltfile  "<< fltfile << endl;
  if ( gvfile.length()>1 && fltfile.length() > 1  ) p4gvf  = new P4gvCor( gvfile, fltfile ) ;


OP's avatar
OP committed
116 117 118
  // ---
  HiStatsInitiator _inia;
  int rc = 0;
OP's avatar
modifs  
OP committed
119 120 121 122
  // Gain correction class
  P4gnuGain * p4g = NULL ;
  try {
    p4g = new P4gnuGain(params.gain_gnu_file_);
123 124
  }
  catch (PException& exc) {
OP's avatar
modifs  
OP committed
125
    cerr << " visiavg.cc catched PException " << exc.Msg() << endl;
126
    cerr << " .... Continuing execution ... " << endl;
OP's avatar
modifs  
OP committed
127 128
    p4g = NULL ;
  }
129 130 131 132 133 134 135 136 137 138 139 140
  P4PhaseCor * p4phc = NULL;
  if (params.phase_corr_param_file_.length() > 0)  {
    try {
      p4phc = new P4PhaseCor(params.phase_corr_param_file_);
    }
    catch (PException& exc) {
      cerr << " visiavg.cc catched PException when reading Phase Correction parameter file: " << params.phase_corr_param_file_ << endl;
      cerr << " Exception message: " << exc.Msg() << endl;
      cerr << " .... Continuing execution ... " << endl;
      p4phc = NULL;
    }
  }
OP's avatar
modifs  
OP committed
141

OP's avatar
OP committed
142 143 144 145 146 147 148 149 150
  try {
    ResourceUsage resu;


    // setting up input visi reader
    vector<string> paths; 
    paths.push_back(params.inpath5_); 
    paths.push_back(params.inpath6_);

151

152
    VisiP4WindowReader wreader(paths, params.inwsz_);
153
    wreader.getReader().setFreqReordering(params.fgreorderfreq_);
OP's avatar
OP committed
154 155
    if (!params.fgserall_ && !params.fgtmsel_) {
      cout << " vreader.SelectSerialNum(Imin="<<Imin<<" ,Imax="<<Imax<<" ,Istep="<<Istep<<")"<<endl;
156
      wreader.getReader().SelectSerialNum(Imin,Imax,Istep);
OP's avatar
OP committed
157 158 159
    }
    else if (params.fgserall_) {
      cout << " vreader.SelectAll() ... " << endl;
160
      wreader.getReader().SelectAll();
OP's avatar
OP committed
161 162 163 164 165 166
    }
    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;
167
      wreader.getReader().SelectTimeFrame(tustart, tuend, Istep);
OP's avatar
OP committed
168
    }
169
    wreader.getReader().setPrintLevel(prtlev);
OP's avatar
OP committed
170

171 172 173
    Imin = wreader.getReader().getSerialFirst();
    Imax =  wreader.getReader().getSerialLast();
    Istep =  wreader.getReader().getSerialStep();
OP's avatar
OP committed
174 175
    cout << "visiavg/Info: processing visibility matrix serial/sequence number range "
	 <<Imin<<" <= seq <= " << Imax << " with step="<<Istep<<endl;
176
    cout << " WindowSize="<<wreader.getWindowSize()<<"  -> TotalNbWindows="<<wreader.getTotalNbWindows()<<endl;
OP's avatar
OP committed
177 178

    bool fgok=true;
OP's avatar
OP committed
179 180
    // vecteur de noms 
    vector <string> ext_names;
OP's avatar
OP committed
181
    // un vecteur avec les temps 
182 183
    TVector< double > timevec(wreader.getTotalNbWindows()/deltaIavg); 
    TVector< double > ravec(wreader.getTotalNbWindows()/deltaIavg); 
OP's avatar
OP committed
184 185
    TMatrix< complex<r_4> > vismtx;
    TMatrix< complex<r_4> > acsum;
OP's avatar
OP committed
186
    TMatrix< r_4 > acsum_sq; // wil sum only the real part ^2
OP's avatar
OP committed
187
    TMatrix< complex<r_4> > cxsum;
OP's avatar
OP committed
188 189 190
    // for sums of real and imag parts 
    TMatrix< r_4 > cxsum_sq_rp;
    TMatrix< r_4 > cxsum_sq_ip;
OP's avatar
OP committed
191 192 193 194 195 196

    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; 
OP's avatar
OP committed
197
    // for
OP's avatar
OP committed
198 199
    //----- 6 H-H cross-cor TimeFrequency maps 
    vector< TArray< complex<r_4> > > vtfm;
OP's avatar
OP committed
200 201 202 203
    //----- 6 H-H cross-cor TimeFrequency maps for the variances of real and imag parts 
    vector< TArray< r_4 > > vtfm_rp_sq;
    vector< TArray< r_4 > > vtfm_ip_sq;

OP's avatar
OP committed
204 205
    //----- 8 auto-corr TimeFrequency maps 
    vector< TArray< r_4 > > vtfmac;
OP's avatar
OP committed
206 207
    //----- 8 auto-corr TimeFrequency variance maps 
    vector< TArray< r_4 > > vtfmac_sq;
OP's avatar
OP committed
208 209 210 211 212 213
    
    //---- for the time-freqency map filling    
    sa_size_t TFMtmidx=0;
    sa_size_t tfmSX, tfmSY;

    while (fgok) {
214
      //reads next visibility matrix window 
215
      fgok = wreader.Shift();
OP's avatar
OP committed
216
      if (!fgok)  break;
217 218 219 220
      
      vismtx = wreader.getAverageVisMtx();
      cfdate = wreader.getAverageTimeTU();
      
OP's avatar
OP committed
221
      // Apply gain g(nu)
OP's avatar
modifs  
OP committed
222 223
      if (p4g != NULL )
	p4g->applyGain(vismtx);
224
      // Apply time dependent (/ temperature) gain G(t) 
OP's avatar
modifs  
OP committed
225 226
      if (p4gvf != NULL) 
	p4gvf->applyGain(vismtx,cfdate);
227 228 229 230 231
      // Apply phase correction
      if (p4phc != NULL) 
	p4phc->applyPhase(vismtx);
      
       if (cnt==0)  {    //resizing matrices for sum of auto-correlations and sum of 6 cross-correlations 
OP's avatar
OP committed
232
	acsum.SetSize(8, vismtx.NCols());
OP's avatar
OP committed
233
	acsum_sq.SetSize(8, vismtx.NCols());
OP's avatar
OP committed
234
	cxsum.SetSize(6, vismtx.NCols());
OP's avatar
OP committed
235 236 237
	cxsum_sq_rp.SetSize(6, vismtx.NCols());
	cxsum_sq_ip.SetSize(6, vismtx.NCols());

238
	tfmSX=wreader.getTotalNbWindows()/deltaIavg;
OP's avatar
OP committed
239 240 241 242
	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) ); 
OP's avatar
OP committed
243 244 245
	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) ); 

OP's avatar
OP committed
246 247 248
	//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) );
OP's avatar
OP committed
249 250 251
	cout << "and the 2x6 for squares of rela & imag 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) );
OP's avatar
OP committed
252 253 254 255 256 257 258 259 260 261

	// recupere le jour de depart @ 0h
	datestart = TimeStamp(cfdate.DaysPart(),0.);
 
      }
      if (I==0) {   // start filling a new time bin 
	dateobs=cfdate;
	if (prtlev>0) 
	  cout<<"visiavg/Info:  dateobs="<<dateobs<<" SecondsPart()="<<dateobs.SecondsPart()<<endl;
	acsum = complex<r_4>(0.,0.);
OP's avatar
OP committed
262
	acsum_sq = 0.;
OP's avatar
OP committed
263
	cxsum = complex<r_4>(0.,0.);
OP's avatar
OP committed
264 265
	cxsum_sq_rp = 0.;
	cxsum_sq_ip = 0.;
OP's avatar
OP committed
266 267 268 269
      }
      
      //   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 
OP's avatar
OP committed
270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286
      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) ;

      }


OP's avatar
OP committed
287 288 289
      I++;    // incrementing DeltaTime counter 
      
      if (I==deltaIavg) {
OP's avatar
OP committed
290
 
OP's avatar
OP committed
291 292 293
	//---- 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));
OP's avatar
OP committed
294
	  TVector<r_4> vacsq = acsum_sq.Row(k);
OP's avatar
OP committed
295 296
	  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];
OP's avatar
OP committed
297
	    TArray< r_4 > & tfmap_sq = vtfmac_sq[k];
OP's avatar
OP committed
298 299
	    for(sa_size_t jy=0; jy<tfmSY; jy++) {  // frequency binning 
	      tfmap(TFMtmidx, jy) = vac( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
OP's avatar
OP committed
300 301
	      tfmap_sq(TFMtmidx, jy) = vacsq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();

OP's avatar
OP committed
302
	    } 
OP's avatar
OP committed
303
	  }
OP's avatar
OP committed
304
	}  //----- end of loop over the 8 AutoCor
OP's avatar
OP committed
305

OP's avatar
OP committed
306
	//---- On s'occupe des 6 cross-correlations  1H-2H ... 3H-4H 
OP's avatar
OP committed
307 308
 	for(size_t k=0; k<KVCXHH.size(); k++)   {   // loop over the 6 Xcor 	  
	  TVector< complex<r_4> > vcx = cxsum.Row(k);
OP's avatar
OP committed
309 310
	  TVector<r_4>  vcxprsq = cxsum_sq_rp.Row(k);
	  TVector<r_4>  vcxpisq = cxsum_sq_ip.Row(k);
OP's avatar
OP committed
311 312
	  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];
OP's avatar
OP committed
313 314
	    TArray< r_4 > & tfmapsqpr = vtfm_rp_sq[k];
	    TArray< r_4 > & tfmapsqpi = vtfm_ip_sq[k];
OP's avatar
OP committed
315 316
	    for(sa_size_t jy=0; jy<tfmSY; jy++) {
	      tfmap(TFMtmidx, jy) = vcx( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
OP's avatar
OP committed
317 318
	      tfmapsqpr(TFMtmidx, jy) = vcxprsq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
	      tfmapsqpi(TFMtmidx, jy) = vcxpisq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
OP's avatar
OP committed
319 320 321
	    } 
	  }  
	}  //----- end of loop over the 6 Xcor 
OP's avatar
OP committed
322 323 324 325
	
	double tdif =  cfdate.TimeDifferenceSeconds(cfdate,dateobs)/2.;
	timevec(TFMtmidx) = dateobs.TimeDifferenceSeconds(dateobs.ShiftSeconds (tdif ),datestart);	// centre du bin 
	ravec(TFMtmidx) = P4Coords::RAFromTimeTU(dateobs);
OP's avatar
OP committed
326 327 328 329 330 331 332 333 334 335 336 337
	TFMtmidx++;
	//  ... done 
	I=0;  cntnt++;
      }
      cnt++;
      if ((cnt>0)&&(cntnt%10==0)&&(cntnt>pcntnt)) {
	cout<<"visiavg/Info: TFM-Map fill cnt="<<cntnt<<" VisMtxCount="<<cnt<<" /Max="<<Imax<<" DateObs="<<dateobs<<endl;
	pcntnt=cntnt;
      }
    }
    cout<<"visiavg/Info: count="<<cnt<<" visimtx read "<<endl;

OP's avatar
OP committed
338
    // --- Sauvegarde cartes temps-frequence en fits 
OP's avatar
OP committed
339 340 341
    //FitsABTWriter * fbtw = NULL;
    FitsInOutFile  * fos = NULL ;

OP's avatar
OP committed
342
    if (fitsoutfile.length()>=1){
OP's avatar
OP committed
343 344 345 346
      cout << " fitsoutfile" <<fitsoutfile<< endl;
      fos = new FitsInOutFile(fitsoutfile, FitsInOutFile::Fits_Create);
    }
    
OP's avatar
OP committed
347

OP's avatar
OP committed
348
    POutPersist potfm(outfile);
OP's avatar
OP committed
349 350
    char bufnam[10];
    int numkey =0;
OP's avatar
OP committed
351 352 353
    // --- renormalizing and saving AutoCorr time-frequency maps 
    cout<<"  visiavg/Info: Saving 8 AutoCorr time-frequency maps to PPF file "<<outfile<<endl;
    const char* tfm_names[8]={"TFM_1H", "TFM_2H", "TFM_3H", "TFM_4H", "TFM_1V", "TFM_2V", "TFM_3V", "TFM_4V"};
OP's avatar
OP committed
354
    const char* tfmsq_names[8]={"VARTFM_1H", "VARTFM_2H", "VARTFM_3H", "VARTFM_4H", "VARTFM_1V", "VARTFM_2V", "VARTFM_3V", "VARTFM_4V"};
OP's avatar
OP committed
355 356 357 358
    for(int k=0; k<8; k++)  {   // loop over the 8 AutoCorr 
      TArray< r_4 > & tfmap = vtfmac[k];
      tfmap *= (r_4)(1./((double)deltaIavg*(double)TFMfbin));
      potfm << PPFNameTag(tfm_names[k]) << tfmap;
OP's avatar
OP committed
359

OP's avatar
OP committed
360 361 362
      if (fos != NULL) {
	ext_names.push_back(tfm_names[k]);
	(*fos)<<  tfmap;
OP's avatar
OP committed
363 364 365 366 367 368
      }
      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) {
OP's avatar
OP committed
369
	ext_names.push_back(tfmsq_names[k]);
OP's avatar
OP committed
370
	(*fos)<<  tfmapsq;
OP's avatar
OP committed
371
      }
OP's avatar
OP committed
372
      
OP's avatar
OP committed
373
    }
OP's avatar
OP committed
374

OP's avatar
OP committed
375 376
    // --- renormalizing and saving H-H Cross-Corr time-frequency maps 
    cout<<"  visiavg/Info: Saving 6 H-H cross-corr time-frequency maps to PPF file "<<outfile<<endl;
OP's avatar
OP committed
377

OP's avatar
OP committed
378
    const char* tfmCC_names[6]={"TFM_1H2H", "TFM_1H3H", "TFM_1H4H", "TFM_2H3H", "TFM_2H4H", "TFM_3H4H"};
OP's avatar
OP committed
379 380 381
    const char* vrtfmCC_names[6]={"RVARTFM_1H2H", "RVARTFM_1H3H", "RVARTFM_1H4H", "RVARTFM_2H3H", "RVARTFM_2H4H", "RVARTFM_3H4H"};
    const char* vitfmCC_names[6]={"IVARTFM_1H2H", "IVARTFM_1H3H", "IVARTFM_1H4H", "IVARTFM_2H3H", "IVARTFM_2H4H", "IVARTFM_3H4H"};

OP's avatar
OP committed
382 383
    for(int k=0; k<6; k++)  {   // loop over the 6 Xcor 
      TArray< complex<r_4> > & tfmap = vtfm[k];
OP's avatar
OP committed
384 385 386

      TArray< r_4 > & tfmap_sqpr = vtfm_rp_sq[k];
      TArray< r_4 > & tfmap_sqpi = vtfm_ip_sq[k];
OP's avatar
OP committed
387
      tfmap *= complex<r_4>((r_4)(1./((double)deltaIavg*(double)TFMfbin)), 0.);
OP's avatar
OP committed
388 389 390 391 392 393 394
      potfm << PPFNameTag(tfmCC_names[k]) << tfmap;
      if (fos != NULL) {
	ext_names.push_back(string(tfmCC_names[k])+"_real");
	(*fos)<<  real(tfmap);
	ext_names.push_back(string(tfmCC_names[k])+"_imag");
	(*fos)<<  imag(tfmap);
      }
OP's avatar
OP committed
395 396 397 398 399 400 401 402 403 404 405
      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;
OP's avatar
OP committed
406 407 408 409 410 411
      if (fos != NULL) {
	ext_names.push_back(vrtfmCC_names[k]);
	(*fos)<<  tfmap_sqpr;
	ext_names.push_back(vitfmCC_names[k]);
	(*fos)<<  tfmap_sqpi;
      }
OP's avatar
OP committed
412

OP's avatar
OP committed
413
    }
OP's avatar
OP committed
414

OP's avatar
OP committed
415
    potfm << PPFNameTag("TimeVec") << timevec ;
OP's avatar
OP committed
416
    potfm << PPFNameTag("RAVec") << ravec ;
OP's avatar
OP committed
417 418
    // --- FIN sauvegarde cartes temps-frequence 
    //    resu.Update();
OP's avatar
OP committed
419 420 421
    P4FreqBand myp4fre;
    TVector <double> avg_freqs( myp4fre.getP4NbFreqChannels()/TFMfbin);
    double frbase =  myp4fre.freqstart_ + myp4fre.getP4FreqResolution()/2. ;
OP's avatar
OP committed
422
    for (int kf=0 ; kf< myp4fre.getP4NbFreqChannels()/TFMfbin ; kf++,frbase += myp4fre.getP4FreqResolution()*TFMfbin )
OP's avatar
OP committed
423
      avg_freqs(kf) = frbase ;
OP's avatar
OP committed
424

OP's avatar
OP committed
425
    potfm << PPFNameTag("FreqVec") << avg_freqs;
OP's avatar
OP committed
426 427 428 429 430 431 432 433 434 435 436
    if (fos != NULL) {
      ext_names.push_back("Frequences");
      (*fos)<< avg_freqs ;
      ext_names.push_back("RAs");
      (*fos)<<  ravec;
      ext_names.push_back("Times");
      (*fos)<<  timevec;
      cout << " number of objs in fits "<< ext_names.size() << endl;
      cout << ext_names << endl;
      delete(fos);
    }
OP's avatar
OP committed
437

OP's avatar
OP committed
438
    cout << resu;   // Update est fait lors du print
OP's avatar
modifs  
OP committed
439
    if (p4g != NULL) delete(p4g);
OP's avatar
OP committed
440
    
OP's avatar
OP committed
441
    
OP's avatar
OP committed
442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462
  }
  catch (PException& exc) {
    cerr << " visiavg.cc catched PException " << exc.Msg() << endl;
    rc = 77;
  }  
  catch (std::exception& sex) {
    cerr << "\n visiavg.cc std::exception :" 
         << (string)typeid(sex).name() << "\n msg= " 
         << sex.what() << endl;
    rc = 78;
  }
  catch (...) {
    cerr << " visiavg.cc catched unknown (...) exception  " << endl; 
    rc = 79; 
  } 

  cout << ">>>> visiavg.cc ------- END ----------- RC=" << rc << endl;
  return rc;

}