visiavg.cc 15 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 43
#include "p4gnugain.h"

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

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

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

OP's avatar
OP committed
69 70 71 72
  string fitsoutfile = params.fitsoutfile_;
  if (fitsoutfile.length()>=1) {
    fitsoutfile = "!"+fitsoutfile ; // adds '!' ?
  }
OP's avatar
OP committed
73 74 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
  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
OP committed
104
  
OP's avatar
OP committed
105 106 107 108 109 110 111 112 113 114 115 116 117 118
  // ---
  HiStatsInitiator _inia;
  int rc = 0;
  try {
    ResourceUsage resu;

    // Gain correction class
    P4gnuGain p4g(params.gain_gnu_file_);

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

119

120
    VisiP4WindowReader wreader(paths, params.inwsz_);
121
    wreader.getReader().setFreqReordering(params.fgreorderfreq_);
OP's avatar
OP committed
122 123
    if (!params.fgserall_ && !params.fgtmsel_) {
      cout << " vreader.SelectSerialNum(Imin="<<Imin<<" ,Imax="<<Imax<<" ,Istep="<<Istep<<")"<<endl;
124
      wreader.getReader().SelectSerialNum(Imin,Imax,Istep);
OP's avatar
OP committed
125 126 127
    }
    else if (params.fgserall_) {
      cout << " vreader.SelectAll() ... " << endl;
128
      wreader.getReader().SelectAll();
OP's avatar
OP committed
129 130 131 132 133 134
    }
    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;
135
      wreader.getReader().SelectTimeFrame(tustart, tuend, Istep);
OP's avatar
OP committed
136
    }
137
    wreader.getReader().setPrintLevel(prtlev);
OP's avatar
OP committed
138

139 140 141
    Imin = wreader.getReader().getSerialFirst();
    Imax =  wreader.getReader().getSerialLast();
    Istep =  wreader.getReader().getSerialStep();
OP's avatar
OP committed
142 143
    cout << "visiavg/Info: processing visibility matrix serial/sequence number range "
	 <<Imin<<" <= seq <= " << Imax << " with step="<<Istep<<endl;
144
    cout << " WindowSize="<<wreader.getWindowSize()<<"  -> TotalNbWindows="<<wreader.getTotalNbWindows()<<endl;
OP's avatar
OP committed
145 146

    bool fgok=true;
OP's avatar
OP committed
147 148
    // vecteur de noms 
    vector <string> ext_names;
OP's avatar
OP committed
149
    // un vecteur avec les temps 
150 151
    TVector< double > timevec(wreader.getTotalNbWindows()/deltaIavg); 
    TVector< double > ravec(wreader.getTotalNbWindows()/deltaIavg); 
OP's avatar
OP committed
152 153
    TMatrix< complex<r_4> > vismtx;
    TMatrix< complex<r_4> > acsum;
OP's avatar
OP committed
154
    TMatrix< r_4 > acsum_sq; // wil sum only the real part ^2
OP's avatar
OP committed
155
    TMatrix< complex<r_4> > cxsum;
OP's avatar
OP committed
156 157 158
    // for sums of real and imag parts 
    TMatrix< r_4 > cxsum_sq_rp;
    TMatrix< r_4 > cxsum_sq_ip;
OP's avatar
OP committed
159 160 161 162 163 164

    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
165
    // for
OP's avatar
OP committed
166 167
    //----- 6 H-H cross-cor TimeFrequency maps 
    vector< TArray< complex<r_4> > > vtfm;
OP's avatar
OP committed
168 169 170 171
    //----- 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
172 173
    //----- 8 auto-corr TimeFrequency maps 
    vector< TArray< r_4 > > vtfmac;
OP's avatar
OP committed
174 175
    //----- 8 auto-corr TimeFrequency variance maps 
    vector< TArray< r_4 > > vtfmac_sq;
OP's avatar
OP committed
176 177 178 179 180 181
    
    //---- for the time-freqency map filling    
    sa_size_t TFMtmidx=0;
    sa_size_t tfmSX, tfmSY;

    while (fgok) {
182
      //reads next visibility matrix window 
183
      fgok = wreader.Shift();
OP's avatar
OP committed
184
      if (!fgok)  break;
185 186 187 188
      
      vismtx = wreader.getAverageVisMtx();
      cfdate = wreader.getAverageTimeTU();
      
OP's avatar
OP committed
189 190
      // Apply gain g(nu)
      p4g.applyGain(vismtx);
OP's avatar
OP committed
191

OP's avatar
OP committed
192

OP's avatar
OP committed
193 194
      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
195
	acsum_sq.SetSize(8, vismtx.NCols());
OP's avatar
OP committed
196
	cxsum.SetSize(6, vismtx.NCols());
OP's avatar
OP committed
197 198 199
	cxsum_sq_rp.SetSize(6, vismtx.NCols());
	cxsum_sq_ip.SetSize(6, vismtx.NCols());

200
	tfmSX=wreader.getTotalNbWindows()/deltaIavg;
OP's avatar
OP committed
201 202 203 204
	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
205 206 207
	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
208 209 210
	//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
211 212 213
	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
214 215 216 217 218 219 220 221 222 223

	// 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
224
	acsum_sq = 0.;
OP's avatar
OP committed
225
	cxsum = complex<r_4>(0.,0.);
OP's avatar
OP committed
226 227
	cxsum_sq_rp = 0.;
	cxsum_sq_ip = 0.;
OP's avatar
OP committed
228 229 230 231
      }
      
      //   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
232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248
      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
249 250 251
      I++;    // incrementing DeltaTime counter 
      
      if (I==deltaIavg) {
OP's avatar
OP committed
252
 
OP's avatar
OP committed
253 254 255
	//---- 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
256
	  TVector<r_4> vacsq = acsum_sq.Row(k);
OP's avatar
OP committed
257 258
	  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
259
	    TArray< r_4 > & tfmap_sq = vtfmac_sq[k];
OP's avatar
OP committed
260 261
	    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
262 263
	      tfmap_sq(TFMtmidx, jy) = vacsq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();

OP's avatar
OP committed
264
	    } 
OP's avatar
OP committed
265
	  }
OP's avatar
OP committed
266
	}  //----- end of loop over the 8 AutoCor
OP's avatar
OP committed
267

OP's avatar
OP committed
268
	//---- On s'occupe des 6 cross-correlations  1H-2H ... 3H-4H 
OP's avatar
OP committed
269 270
 	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
271 272
	  TVector<r_4>  vcxprsq = cxsum_sq_rp.Row(k);
	  TVector<r_4>  vcxpisq = cxsum_sq_ip.Row(k);
OP's avatar
OP committed
273 274
	  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
275 276
	    TArray< r_4 > & tfmapsqpr = vtfm_rp_sq[k];
	    TArray< r_4 > & tfmapsqpi = vtfm_ip_sq[k];
OP's avatar
OP committed
277 278
	    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
279 280
	      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
281 282 283
	    } 
	  }  
	}  //----- end of loop over the 6 Xcor 
OP's avatar
OP committed
284 285 286 287
	
	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
288 289 290 291 292 293 294 295 296 297 298 299
	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
300
    // --- Sauvegarde cartes temps-frequence en fits 
OP's avatar
OP committed
301 302 303
    //FitsABTWriter * fbtw = NULL;
    FitsInOutFile  * fos = NULL ;

OP's avatar
OP committed
304
    if (fitsoutfile.length()>=1){
OP's avatar
OP committed
305 306 307 308
      cout << " fitsoutfile" <<fitsoutfile<< endl;
      fos = new FitsInOutFile(fitsoutfile, FitsInOutFile::Fits_Create);
    }
    
OP's avatar
OP committed
309

OP's avatar
OP committed
310
    POutPersist potfm(outfile);
OP's avatar
OP committed
311 312
    char bufnam[10];
    int numkey =0;
OP's avatar
OP committed
313 314 315
    // --- 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
316
    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
317 318 319 320
    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
321

OP's avatar
OP committed
322 323 324
      if (fos != NULL) {
	ext_names.push_back(tfm_names[k]);
	(*fos)<<  tfmap;
OP's avatar
OP committed
325 326 327 328 329 330
      }
      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
331
	ext_names.push_back(tfmsq_names[k]);
OP's avatar
OP committed
332
	(*fos)<<  tfmapsq;
OP's avatar
OP committed
333
      }
OP's avatar
OP committed
334
      
OP's avatar
OP committed
335
    }
OP's avatar
OP committed
336

OP's avatar
OP committed
337 338
    // --- 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
339

OP's avatar
OP committed
340
    const char* tfmCC_names[6]={"TFM_1H2H", "TFM_1H3H", "TFM_1H4H", "TFM_2H3H", "TFM_2H4H", "TFM_3H4H"};
OP's avatar
OP committed
341 342 343
    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
344 345
    for(int k=0; k<6; k++)  {   // loop over the 6 Xcor 
      TArray< complex<r_4> > & tfmap = vtfm[k];
OP's avatar
OP committed
346 347 348

      TArray< r_4 > & tfmap_sqpr = vtfm_rp_sq[k];
      TArray< r_4 > & tfmap_sqpi = vtfm_ip_sq[k];
OP's avatar
OP committed
349
      tfmap *= complex<r_4>((r_4)(1./((double)deltaIavg*(double)TFMfbin)), 0.);
OP's avatar
OP committed
350 351 352 353 354 355 356
      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
357 358 359 360 361 362 363 364 365 366 367
      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
368 369 370 371 372 373
      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
374

OP's avatar
OP committed
375
    }
OP's avatar
OP committed
376

OP's avatar
OP committed
377
    potfm << PPFNameTag("TimeVec") << timevec ;
OP's avatar
OP committed
378
    potfm << PPFNameTag("RAVec") << ravec ;
OP's avatar
OP committed
379 380
    // --- FIN sauvegarde cartes temps-frequence 
    //    resu.Update();
OP's avatar
OP committed
381 382 383
    P4FreqBand myp4fre;
    TVector <double> avg_freqs( myp4fre.getP4NbFreqChannels()/TFMfbin);
    double frbase =  myp4fre.freqstart_ + myp4fre.getP4FreqResolution()/2. ;
OP's avatar
OP committed
384
    for (int kf=0 ; kf< myp4fre.getP4NbFreqChannels()/TFMfbin ; kf++,frbase += myp4fre.getP4FreqResolution()*TFMfbin )
OP's avatar
OP committed
385
      avg_freqs(kf) = frbase ;
OP's avatar
OP committed
386

OP's avatar
OP committed
387
    potfm << PPFNameTag("FreqVec") << avg_freqs;
OP's avatar
OP committed
388 389 390 391 392 393 394 395 396 397 398
    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
399

OP's avatar
OP committed
400
    cout << resu;   // Update est fait lors du print
OP's avatar
OP committed
401
    
OP's avatar
OP committed
402
    
OP's avatar
OP committed
403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423
  }
  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;

}