visi2ntac.cc 7.95 KB
Newer Older
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  2008-2015 

   visi2ntac: programme de lecture des fichiers matrices de 
   visibilites de PAON4, remplissage et ecriture d'un 
   DataTable avec la puissance des 8 voies en fonction du temps 
   R. Ansari, J.E. Campagne, C. Magneville   -  LAL/Irfu
   ---------------------------------------------------------- */

// 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 "visip4reader.h"
41 42
#include "p4autils.h"
#include "p4gnugain.h"
43 44

//--------------------------- Fonctions de ce fichier   ------------------- 
45 46
int Usage(void);
int Usage(void)
47
{
48 49 50 51 52
  cout << " --- visi2ntac.cc : Read PPF files produced by mfacq and fills an NTuple \n" << endl;
  cout << " Usage: visi2ntac [-arguments] [-cg] [freq1=1393] [freq2=1420.4] \n"
       << "    -cg : apply gain corrections \n"<<endl;
  P4AnaParams::UsageOptions();
  cout<< endl;
53 54 55 56 57
  return 1;
}

//----------------------------------------------------
//----------------------------------------------------
58
int main(int narg, const char* arg[])
59
{
60 61 62 63 64 65 66 67 68 69 70 71 72 73 74
  // --- Decoding parameters 
  if( (narg<2) || ((narg>1)&&(strcmp(arg[1],"-h")==0) ) )  return Usage();
  P4AnaParams params;
  params.DecodeArgs(narg, arg);
  string outfile = params.outfile_;
  if (outfile.length()<1)  outfile = "visdt.ppf";
  int Imin = params.Imin_, Imax = params.Imax_, Istep = params.Istep_; 
  bool fgcorgain = false;
  double freq1 = 1393.;
  double freq2 = 1420.4;
  if (params.lastargs_.size()>0) {
    size_t off=0;
    if (params.lastargs_[0] == "-cg")  { fgcorgain=true;   off=1; }
    if (params.lastargs_.size() > off)  freq1 = atof(params.lastargs_[off].c_str());
    if (params.lastargs_.size() > off+1)  freq2 = atof(params.lastargs_[off+1].c_str());
75
  }
76 77 78 79 80 81 82 83 84 85 86 87 88 89 90
  P4FreqBand fb1(freq1,0.);
  sa_size_t JF1 =  fb1.getFirstFreqChannel();
  P4FreqBand fb2(freq1,0.);
  sa_size_t JF2 =  fb1.getFirstFreqChannel();

  params.Print(cout);
  cout <<"visi2ntac/Info: Path BAO5:"<<params.inpath5_<<" BAO6:"<<params.inpath6_<<"\n"
       <<"fgreorderfreq="<<params.fgreorderfreq_<<"\n"
       <<"Imin,max,step="<<Imin<<","<<Imax<<","<<Istep
       <<"outfile="<<outfile<<endl;
  cout << (fgcorgain ? " APPLY GainCorrection " : " NO GainCorr") << " Freq1="<<freq1<<" MHz -> JF1="<<JF1
       <<" Freq1="<<freq1<<" MHz -> JF1="<<JF1<<endl;

  P4gnuGain * p4gp=NULL;
  if (fgcorgain)  p4gp = new  P4gnuGain(params.gain_gnu_file_);
91
 
92 93 94 95 96 97
  HiStatsInitiator _inia;
  //   TArrayInitiator  _inia;

  int rc = 0;
  try {
    ResourceUsage resu;
98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113

    const char* dtnames[51]={ "k","day","time",
			      "ac1h","ac2h","ac3h","ac4h","ac1v","ac2v","ac3v","ac4v",
			      "ac1hB","ac2hB","ac3hB","ac4hB","ac1vB","ac2vB","ac3vB","ac4vB",
			      "cx1h2h","cx1h3h","cx1h4h","cx2h3h","cx2h4h","cx3h4h",
			      "cx1h2hB","cx1h3hB","cx1h4hB","cx2h3hB","cx2h4hB","cx3h4hB",
			      "cx1v2v","cx1v3v","cx1v4v","cx2v3v","cx2v4v","cx3v4v",
			      "cx1v2vB","cx1v3vB","cx1v4vB","cx2v3vB","cx2v4vB","cx3v4vB",
			      "cx1h4v","cx2h4v","cx3h4v","cx4h4v",
			      "cx1h4vB","cx2h4vB","cx3h4vB","cx4h4vB" };


    DataTable dt(384);
    dt.AddIntegerColumn(dtnames[0]);
    dt.AddIntegerColumn(dtnames[1]);
    dt.AddFloatColumn(dtnames[2]);
114
    for(int kk=0; kk<16; kk++) 
115 116
      dt.AddFloatColumn(dtnames[3+kk]);
    for(int kk=0; kk<24; kk++) 
117
      dt.AddComplexColumn(dtnames[3+16+kk]);
118
    for(int kk=0; kk<8; kk++) 
119
      dt.AddComplexColumn(dtnames[3+16+24+kk]);
120 121 122 123 124 125 126 127 128 129 130

    P4AVisiNumEncoder p4venc;
    std::vector<sa_size_t> kac = p4venc.getAllAutoCor();
    std::vector<sa_size_t> kcxh = p4venc.getAllHCrossCor();
    std::vector<sa_size_t> kcxv = p4venc.getAllVCrossCor();
    std::vector<sa_size_t> kcxh4v(4);
    kcxh4v[0] = p4venc.Convert2VisiNumber("1H4V");
    kcxh4v[1] = p4venc.Convert2VisiNumber("2H4V");
    kcxh4v[2] = p4venc.Convert2VisiNumber("3H4V");
    kcxh4v[3] = p4venc.Convert2VisiNumber("4H4V");
    
131
    vector<string> paths; 
132 133
    paths.push_back(params.inpath5_); 
    paths.push_back(params.inpath6_);
134

135 136
    VisiP4ReaderBase * reader = VisiP4ReaderBase::getReader(paths);
    VisiP4ReaderBase & vreader = (*reader);
137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157
    vreader.setFreqReordering(params.fgreorderfreq_);
    if (!params.fgserall_ && !params.fgtmsel_) {
      cout << " vreader.SelectSerialNum(Imin="<<Imin<<" ,Imax="<<Imax<<" ,Istep="<<Istep<<")"<<endl;
      vreader.SelectSerialNum(Imin,Imax,Istep);
    }
    else if (params.fgserall_) {
      cout << " vreader.SelectAll() ... " << endl;
      vreader.SelectAll();
    }
    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;
      vreader.SelectTimeFrame(tustart, tuend, Istep);
    }
    vreader.setPrintLevel(params.prtlev_);

    Imin = vreader.getSerialFirst(); Imax =  vreader.getSerialLast();   Istep =  vreader.getSerialStep();
    cout << "visi2dtacx/Info: processing visibility matrix serial/sequence number range "
	 <<Imin<<" <= seq <= " << Imax << " with step="<<Istep<<endl;
158
   
159 160
    bool fgok=true;
    TimeStamp dateobs;
161 162 163 164 165 166
    TimeStamp dateorg(2017,9,1,12,0,0.);  // Date origine 1 sep 2017
    int cnt=0;
    
    // Getting en empty row_ptr
    DataTableRowPtr rowp = dt.EmptyRowPtr();
    TMatrix< complex<r_4> > vismtx;
167
    double mttag;
168

169 170
    while (fgok) {
      fgok=vreader.ReadNext(vismtx, dateobs, mttag);
171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203
      if (! fgok) break;
      if (fgcorgain)  p4gp->applyGain(vismtx);
      // remplissage datatable 
      dt.NextRowPtr(rowp);
      rowp(0)=(int_4)cnt;   
      rowp(1)=(int_4)(dateobs.DaysPart()-dateorg.DaysPart());   
      rowp(2)=(r_4)dateobs.SecondsPart();
      size_t off = 3;
      for(size_t k=0; k<kac.size(); k++) {
	rowp(off+k) = (float)(vismtx(kac[k], JF1).real());
	rowp(off+k+kac.size()) = (float)(vismtx(kac[k], JF2).real());
      }
      
      off += (2*kac.size());
      for(size_t k=0; k<kcxh.size(); k++)  {
	rowp(off+k) = complex<float>(vismtx(kcxh[k], JF1));
	rowp(off+k+kcxh.size()) = complex<float>(vismtx(kcxh[k], JF2));
      }
      off += (2*kcxh.size());
      for(size_t k=0; k<kcxv.size(); k++)  {
	rowp(off+k) = complex<float>(vismtx(kcxv[k], JF1));
	rowp(off+k+kcxv.size()) = complex<float>(vismtx(kcxv[k], JF2));
      }
      off += (2*kcxv.size());
      for(size_t k=0; k<kcxh4v.size(); k++)  {
	rowp(off+k) = complex<float>(vismtx(kcxh4v[k], JF1));
	rowp(off+k+kcxh4v.size()) = complex<float>(vismtx(kcxh4v[k], JF2));
      }
      
      cnt++;
      if (cnt%50==0) 
	cout << " visi2ntac/Info: fill-count="<<cnt<<" DateObs="<<dateobs<<"  /Max="<<(Imax-Imin+1)/Istep<<endl;
      
204 205
    }
    cout << " visi2ntac/Info: count="<<cnt<<" visimtx read "<<endl;
206 207
    cout << dt;
    cout<<"  rdvisip4/Info: Saving dt DataTable to PPF file "<<outfile<<endl;
208
    POutPersist po(outfile);
209
    po<<dt; 
210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233
    //    resu.Update();
    cout << resu;   // Update est fait lors du print
  }
  catch (PException& exc) {
    cerr << " visi2ntac.cc catched PException " << exc.Msg() << endl;
    rc = 77;
  }  
  catch (std::exception& sex) {
    cerr << "\n visi2ntac.cc std::exception :" 
         << (string)typeid(sex).name() << "\n msg= " 
         << sex.what() << endl;
    rc = 78;
  }
  catch (...) {
    cerr << " visi2ntac.cc catched unknown (...) exception  " << endl; 
    rc = 79; 
  } 

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

}