visiavg.cc 15.1 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 121

    VisiP4WindowReader wreader(paths);
    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) {
OP's avatar
OP committed
182
      //reads next visimtx
183 184
      //      fgok=vreader.ReadNext(vismtx, cfdate, mttag);
      fgok = wreader.Shift();
OP's avatar
OP committed
185
      if (!fgok)  break;
186 187 188 189
      
      vismtx = wreader.getAverageVisMtx();
      cfdate = wreader.getAverageTimeTU();
      
OP's avatar
OP committed
190 191
      // Apply gain g(nu)
      p4g.applyGain(vismtx);
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 252 253 254
      I++;    // incrementing DeltaTime counter 
      
      if (I==deltaIavg) {
	//---- 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
255
	  TVector<r_4> vacsq = acsum_sq.Row(k);
OP's avatar
OP committed
256 257
	  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
258
	    TArray< r_4 > & tfmap_sq = vtfmac_sq[k];
OP's avatar
OP committed
259 260
	    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
261 262
	      tfmap_sq(TFMtmidx, jy) = vacsq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();

OP's avatar
OP committed
263
	    } 
OP's avatar
OP committed
264 265
	  }
  
OP's avatar
OP committed
266
	}  //----- end of loop over the 8 AutoCor
OP's avatar
OP committed
267
	//---- On s'occupe des 6 cross-correlations  1H-2H ... 3H-4H 
OP's avatar
OP committed
268 269
 	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
270 271
	  TVector<r_4>  vcxprsq = cxsum_sq_rp.Row(k);
	  TVector<r_4>  vcxpisq = cxsum_sq_ip.Row(k);
OP's avatar
OP committed
272 273
	  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
274 275
	    TArray< r_4 > & tfmapsqpr = vtfm_rp_sq[k];
	    TArray< r_4 > & tfmapsqpi = vtfm_ip_sq[k];
OP's avatar
OP committed
276 277
	    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
278 279
	      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
280 281 282
	    } 
	  }  
	}  //----- end of loop over the 6 Xcor 
OP's avatar
OP committed
283 284 285 286
	
	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
287 288 289 290 291 292 293 294 295 296 297 298
	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
299
    // --- Sauvegarde cartes temps-frequence en fits 
OP's avatar
OP committed
300 301 302 303 304 305 306 307 308 309 310 311
    //FitsABTWriter * fbtw = NULL;
    FitsInOutFile  * fos = NULL ;
    if (fitsoutfile.length()>=1){
      //cout<<">>>>> Creating a new fits file with FitsABTWriter"<<endl;
      //fbtw = new FitsABTWriter(fitsoutfile,BINARY_TBL,3);
      //fbtw->SetExtName("Test fits table 1");
      //i1 = fbtw.AddCol("Time",NULL,"unit1",TDOUBLE);

      cout << " fitsoutfile" <<fitsoutfile<< endl;
      fos = new FitsInOutFile(fitsoutfile, FitsInOutFile::Fits_Create);
    }
    
OP's avatar
OP committed
312

OP's avatar
OP committed
313
    POutPersist potfm(outfile);
OP's avatar
OP committed
314 315
    char bufnam[10];
    int numkey =0;
OP's avatar
OP committed
316 317 318
    // --- 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
319
    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
320 321 322
    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));
OP's avatar
OP committed
323 324 325
      TArray< r_4 > & tfmapsq = vtfmac_sq[k];
      tfmapsq *= (r_4)(1./((double)deltaIavg*(double)TFMfbin));
      tfmapsq = tfmapsq - tfmap.MulElt(tfmap,tfmap) ; 
OP's avatar
OP committed
326
      potfm << PPFNameTag(tfm_names[k]) << tfmap;
OP's avatar
OP committed
327
      potfm << PPFNameTag(tfmsq_names[k]) << tfmapsq;
OP's avatar
OP committed
328 329 330 331
      if (fos != NULL) {
	ext_names.push_back(tfm_names[k]);
	(*fos)<<  tfmap;
	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 336 337
    }
    // --- 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
338

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

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

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

OP's avatar
OP committed
383
    potfm << PPFNameTag("FreqVec") << avg_freqs;
OP's avatar
OP committed
384 385 386 387 388 389 390 391 392 393 394
    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
395

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

}