visiavg.cc 15.6 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"
OP's avatar
OP committed
44

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

OP's avatar
OP committed
48 49 50 51 52 53 54 55 56 57 58 59 60 61 62
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
63
  FitsIOServerInit();
OP's avatar
OP committed
64 65
  P4AnaParams params;
  params.DecodeArgs(narg, arg);
OP's avatar
OP committed
66

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

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

OP's avatar
OP committed
75 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
  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
106 107 108 109 110 111 112 113 114

  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
115 116 117
  // ---
  HiStatsInitiator _inia;
  int rc = 0;
OP's avatar
modifs  
OP committed
118 119 120 121 122 123 124 125 126 127
  // Gain correction class
  P4gnuGain * p4g = NULL ;
  try {
    p4g = new P4gnuGain(params.gain_gnu_file_);
  }   catch (PException& exc) {
    cerr << " visiavg.cc catched PException " << exc.Msg() << endl;
    cerr << " on continue quand meme "<<endl ;
    p4g = NULL ;
  }

OP's avatar
OP committed
128 129 130 131 132 133 134 135 136
  try {
    ResourceUsage resu;


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

137

138
    VisiP4WindowReader wreader(paths, params.inwsz_);
139
    wreader.getReader().setFreqReordering(params.fgreorderfreq_);
OP's avatar
OP committed
140 141
    if (!params.fgserall_ && !params.fgtmsel_) {
      cout << " vreader.SelectSerialNum(Imin="<<Imin<<" ,Imax="<<Imax<<" ,Istep="<<Istep<<")"<<endl;
142
      wreader.getReader().SelectSerialNum(Imin,Imax,Istep);
OP's avatar
OP committed
143 144 145
    }
    else if (params.fgserall_) {
      cout << " vreader.SelectAll() ... " << endl;
146
      wreader.getReader().SelectAll();
OP's avatar
OP committed
147 148 149 150 151 152
    }
    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;
153
      wreader.getReader().SelectTimeFrame(tustart, tuend, Istep);
OP's avatar
OP committed
154
    }
155
    wreader.getReader().setPrintLevel(prtlev);
OP's avatar
OP committed
156

157 158 159
    Imin = wreader.getReader().getSerialFirst();
    Imax =  wreader.getReader().getSerialLast();
    Istep =  wreader.getReader().getSerialStep();
OP's avatar
OP committed
160 161
    cout << "visiavg/Info: processing visibility matrix serial/sequence number range "
	 <<Imin<<" <= seq <= " << Imax << " with step="<<Istep<<endl;
162
    cout << " WindowSize="<<wreader.getWindowSize()<<"  -> TotalNbWindows="<<wreader.getTotalNbWindows()<<endl;
OP's avatar
OP committed
163 164

    bool fgok=true;
OP's avatar
OP committed
165 166
    // vecteur de noms 
    vector <string> ext_names;
OP's avatar
OP committed
167
    // un vecteur avec les temps 
168 169
    TVector< double > timevec(wreader.getTotalNbWindows()/deltaIavg); 
    TVector< double > ravec(wreader.getTotalNbWindows()/deltaIavg); 
OP's avatar
OP committed
170 171
    TMatrix< complex<r_4> > vismtx;
    TMatrix< complex<r_4> > acsum;
OP's avatar
OP committed
172
    TMatrix< r_4 > acsum_sq; // wil sum only the real part ^2
OP's avatar
OP committed
173
    TMatrix< complex<r_4> > cxsum;
OP's avatar
OP committed
174 175 176
    // for sums of real and imag parts 
    TMatrix< r_4 > cxsum_sq_rp;
    TMatrix< r_4 > cxsum_sq_ip;
OP's avatar
OP committed
177 178 179 180 181 182

    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
183
    // for
OP's avatar
OP committed
184 185
    //----- 6 H-H cross-cor TimeFrequency maps 
    vector< TArray< complex<r_4> > > vtfm;
OP's avatar
OP committed
186 187 188 189
    //----- 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
190 191
    //----- 8 auto-corr TimeFrequency maps 
    vector< TArray< r_4 > > vtfmac;
OP's avatar
OP committed
192 193
    //----- 8 auto-corr TimeFrequency variance maps 
    vector< TArray< r_4 > > vtfmac_sq;
OP's avatar
OP committed
194 195 196 197 198 199
    
    //---- for the time-freqency map filling    
    sa_size_t TFMtmidx=0;
    sa_size_t tfmSX, tfmSY;

    while (fgok) {
200
      //reads next visibility matrix window 
201
      fgok = wreader.Shift();
OP's avatar
OP committed
202
      if (!fgok)  break;
203 204 205 206
      
      vismtx = wreader.getAverageVisMtx();
      cfdate = wreader.getAverageTimeTU();
      
OP's avatar
OP committed
207
      // Apply gain g(nu)
OP's avatar
modifs  
OP committed
208 209
      if (p4g != NULL )
	p4g->applyGain(vismtx);
OP's avatar
OP committed
210

OP's avatar
modifs  
OP committed
211 212
      if (p4gvf != NULL) 
	p4gvf->applyGain(vismtx,cfdate);
OP's avatar
OP committed
213

OP's avatar
OP committed
214 215
      if (cnt==0)  {    //resizing matrices for sum of auto-correlations and sum of 6 cross-correlations 
	acsum.SetSize(8, vismtx.NCols());
OP's avatar
OP committed
216
	acsum_sq.SetSize(8, vismtx.NCols());
OP's avatar
OP committed
217
	cxsum.SetSize(6, vismtx.NCols());
OP's avatar
OP committed
218 219 220
	cxsum_sq_rp.SetSize(6, vismtx.NCols());
	cxsum_sq_ip.SetSize(6, vismtx.NCols());

221
	tfmSX=wreader.getTotalNbWindows()/deltaIavg;
OP's avatar
OP committed
222 223 224 225
	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
226 227 228
	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
229 230 231
	//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
232 233 234
	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
235 236 237 238 239 240 241 242 243 244

	// 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
245
	acsum_sq = 0.;
OP's avatar
OP committed
246
	cxsum = complex<r_4>(0.,0.);
OP's avatar
OP committed
247 248
	cxsum_sq_rp = 0.;
	cxsum_sq_ip = 0.;
OP's avatar
OP committed
249 250 251 252
      }
      
      //   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
253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269
      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
270 271 272
      I++;    // incrementing DeltaTime counter 
      
      if (I==deltaIavg) {
OP's avatar
OP committed
273
 
OP's avatar
OP committed
274 275 276
	//---- 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
277
	  TVector<r_4> vacsq = acsum_sq.Row(k);
OP's avatar
OP committed
278 279
	  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
280
	    TArray< r_4 > & tfmap_sq = vtfmac_sq[k];
OP's avatar
OP committed
281 282
	    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
283 284
	      tfmap_sq(TFMtmidx, jy) = vacsq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();

OP's avatar
OP committed
285
	    } 
OP's avatar
OP committed
286
	  }
OP's avatar
OP committed
287
	}  //----- end of loop over the 8 AutoCor
OP's avatar
OP committed
288

OP's avatar
OP committed
289
	//---- On s'occupe des 6 cross-correlations  1H-2H ... 3H-4H 
OP's avatar
OP committed
290 291
 	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
292 293
	  TVector<r_4>  vcxprsq = cxsum_sq_rp.Row(k);
	  TVector<r_4>  vcxpisq = cxsum_sq_ip.Row(k);
OP's avatar
OP committed
294 295
	  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
296 297
	    TArray< r_4 > & tfmapsqpr = vtfm_rp_sq[k];
	    TArray< r_4 > & tfmapsqpi = vtfm_ip_sq[k];
OP's avatar
OP committed
298 299
	    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
300 301
	      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
302 303 304
	    } 
	  }  
	}  //----- end of loop over the 6 Xcor 
OP's avatar
OP committed
305 306 307 308
	
	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
309 310 311 312 313 314 315 316 317 318 319 320
	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
321
    // --- Sauvegarde cartes temps-frequence en fits 
OP's avatar
OP committed
322 323 324
    //FitsABTWriter * fbtw = NULL;
    FitsInOutFile  * fos = NULL ;

OP's avatar
OP committed
325
    if (fitsoutfile.length()>=1){
OP's avatar
OP committed
326 327 328 329
      cout << " fitsoutfile" <<fitsoutfile<< endl;
      fos = new FitsInOutFile(fitsoutfile, FitsInOutFile::Fits_Create);
    }
    
OP's avatar
OP committed
330

OP's avatar
OP committed
331
    POutPersist potfm(outfile);
OP's avatar
OP committed
332 333
    char bufnam[10];
    int numkey =0;
OP's avatar
OP committed
334 335 336
    // --- 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
337
    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
338 339 340 341
    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
342

OP's avatar
OP committed
343 344 345
      if (fos != NULL) {
	ext_names.push_back(tfm_names[k]);
	(*fos)<<  tfmap;
OP's avatar
OP committed
346 347 348 349 350 351
      }
      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
352
	ext_names.push_back(tfmsq_names[k]);
OP's avatar
OP committed
353
	(*fos)<<  tfmapsq;
OP's avatar
OP committed
354
      }
OP's avatar
OP committed
355
      
OP's avatar
OP committed
356
    }
OP's avatar
OP committed
357

OP's avatar
OP committed
358 359
    // --- 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
360

OP's avatar
OP committed
361
    const char* tfmCC_names[6]={"TFM_1H2H", "TFM_1H3H", "TFM_1H4H", "TFM_2H3H", "TFM_2H4H", "TFM_3H4H"};
OP's avatar
OP committed
362 363 364
    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
365 366
    for(int k=0; k<6; k++)  {   // loop over the 6 Xcor 
      TArray< complex<r_4> > & tfmap = vtfm[k];
OP's avatar
OP committed
367 368 369

      TArray< r_4 > & tfmap_sqpr = vtfm_rp_sq[k];
      TArray< r_4 > & tfmap_sqpi = vtfm_ip_sq[k];
OP's avatar
OP committed
370
      tfmap *= complex<r_4>((r_4)(1./((double)deltaIavg*(double)TFMfbin)), 0.);
OP's avatar
OP committed
371 372 373 374 375 376 377
      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
378 379 380 381 382 383 384 385 386 387 388
      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
389 390 391 392 393 394
      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
395

OP's avatar
OP committed
396
    }
OP's avatar
OP committed
397

OP's avatar
OP committed
398
    potfm << PPFNameTag("TimeVec") << timevec ;
OP's avatar
OP committed
399
    potfm << PPFNameTag("RAVec") << ravec ;
OP's avatar
OP committed
400 401
    // --- FIN sauvegarde cartes temps-frequence 
    //    resu.Update();
OP's avatar
OP committed
402 403 404
    P4FreqBand myp4fre;
    TVector <double> avg_freqs( myp4fre.getP4NbFreqChannels()/TFMfbin);
    double frbase =  myp4fre.freqstart_ + myp4fre.getP4FreqResolution()/2. ;
OP's avatar
OP committed
405
    for (int kf=0 ; kf< myp4fre.getP4NbFreqChannels()/TFMfbin ; kf++,frbase += myp4fre.getP4FreqResolution()*TFMfbin )
OP's avatar
OP committed
406
      avg_freqs(kf) = frbase ;
OP's avatar
OP committed
407

OP's avatar
OP committed
408
    potfm << PPFNameTag("FreqVec") << avg_freqs;
OP's avatar
OP committed
409 410 411 412 413 414 415 416 417 418 419
    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
420

OP's avatar
OP committed
421
    cout << resu;   // Update est fait lors du print
OP's avatar
modifs  
OP committed
422
    if (p4g != NULL) delete(p4g);
OP's avatar
OP committed
423
    
OP's avatar
OP committed
424
    
OP's avatar
OP committed
425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445
  }
  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;

}