visi2dtacx.cc 13.9 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 41 42 43 44 45 46 47 48 49 50
// Utilisation de SOPHYA pour faciliter les tests ...
#include "sopnamsp.h"
#include "machdefs.h"

/* ---------------------------------------------------------- 
   Projet BAORadio/PAON4 - (C) LAL/IRFU  2008-2015 

   visi2dtacx: programme de lecture des fichiers matrices de 
   visibilites de PAON4, remplissage et ecriture d'un 
   DataTable avec les cross-correlations 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"

//--------------------------- Fonctions de ce fichier   ------------------- 
int Usage(bool fgshort=true);
// int DecodeArgs(int narg, char* arg[]);

/* --Fonction-- */
int Usage(bool fgea)
{
  cout << " --- visi2dtacx.cc : Read PPF files produced by mfacq time-frequency\n" << endl;
  if (fgea) cout << " visi2dtacx: Missing or wrong argument ! " << endl;
51 52
  cout << " Usage: visi2dtacx [-tfm Nf] InPathBAO5 InPathBAO6 Imin,Imax GainPPFFile OutPPFFile Freq1 Freq2 DeltaFreq DeltaIAvg [PrtLev=0] \n"
       << "  -tfm Nf : create time-frequncy maps of cross-correlation binning Nf frequencies \n" 
53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
       << "  InPathBAO5/6: path for input vismtx_J_iii.ppf on bao5/6 \n"
       << "  Imin,Imax : read files for Imin<=iii<=Imax  iii+=step \n"
       << "  GainPPFFile: Input gains PPF file name \n"
       << "  OutPPFFile: Output PPF file name \n"
       << "  Freq1 Freq2: the two frequencies for which cross-correlations are computed [=1375 1410 MHz] \n"
       << "  DeltaFreq: frequency band (in MHz) [=1 MHz] \n"
       << "  DeltaIAvg: compute average power every DeltaI input vismtx files, def=10 \n"
       << "    Note: bande_freq=[freq1,2-deltafreq , freq1,2+deltafreq] \n"<< endl;
  /*
  if (fgshort) {
    cout << " rdvisip4 -h for detailed instructions" << endl;
    return 1;
  }
  */
  return 1;
}

//----------------------------------------------------
//----------------------------------------------------
int main(int narg, char* arg[])
{
  if ((narg>1)&&(strcmp(arg[1],"-h")==0))  return Usage(false);
  if (narg<10) return Usage(true);
76 77 78 79 80 81 82 83 84 85 86 87
  int offa=0;
  bool FgTFM=false;   // true -> create time-frequency maps of cross-correlations 
  sa_size_t TFMfbin=8;
  if (strcmp(arg[1],"-tfm")==0) {
    if (narg<12) return Usage(true);
    FgTFM=true;
    TFMfbin=atoi(arg[2]);
    if (TFMfbin<1) TFMfbin=8;
    offa=2;
  }
  string path5 = arg[offa+1];
  string path6 = arg[offa+2];
88
  int Imin,Imax,Istep=1; 
89
  sscanf(arg[offa+3],"%d,%d",&Imin,&Imax);
90
  Istep=1;
91 92
  string gainfile=arg[offa+4];
  string outfile=arg[offa+5];
93 94 95 96
  // frequency range definition 
  double deltafreq=1.;
  double freq1=1375.;
  double freq2=1410.;
97 98 99 100
  int nargo=narg-offa;
  if ((nargo>7)&&(strcmp(arg[offa+6],"-")!=0))  freq1=atof(arg[offa+6]);
  if ((nargo>8)&&(strcmp(arg[offa+7],"-")!=0))  freq2=atof(arg[offa+7]);
  if ((nargo>6)&&(strcmp(arg[offa+8],"-")!=0))  deltafreq=atof(arg[offa+8]);
101 102
  // time averaging interval definition, in term of visibility files
  int deltaIavg=10;
103
  if ((nargo>9)&&(strcmp(arg[offa+9],"-")!=0))  deltaIavg=atoi(arg[offa+9]);
104 105
  if (deltaIavg<1)  deltaIavg=1;
  int prtlev=1;
106
  if (nargo>10) prtlev=atoi(arg[offa+10]);
107 108 109 110 111 112 113

  cout << " visi2dtacx/Info: Path BAO5:"<<path5<<" BAO6:"<<path6<<"\n"
       << " Imin,max,step="<<Imin<<","<<Imax<<","<<Istep<<"  OutFile:"<<outfile<<"\n"
       << " freq1="<<freq1<<" freq2="<<freq2<<" deltafreq="<<deltafreq<<" MHz\n"
    //       << JFmin<<" <=NumFreq<= "<<JFmax<<"  "<<JFminB<<" <=NumFreqB<= "<<JFmaxB
       <<"  DeltaIAvg="<<deltaIavg<<" PrtLev="<<prtlev<<endl;

114 115 116 117 118
  string TFMoutfile = "tfm_"+outfile;
  if (FgTFM) {
    cout << " visi2dtacx/Info: Time-Frequency maps of Xcor will be created and saved to file "
	 <<TFMoutfile<<endl;
  }
119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 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 204 205 206 207 208 209 210 211

  //---- calcul des index de bandes de frequences 
  // pas en frequence du firmware FFT 
  double deltanufft=250./4096;    // 250 MHz en 4096 frequences 
  sa_size_t DJF=deltafreq/deltanufft;
  double freqstart=1250.;   // Bande de 1250-1500 MHz 
  sa_size_t JFmin=(freq1-1250.-deltafreq)/deltanufft;
  sa_size_t JFmax=(freq1-1250.+deltafreq)/deltanufft;
  sa_size_t JFminB=(freq2-1250.-deltafreq)/deltanufft;
  sa_size_t JFmaxB=(freq2-1250.+deltafreq)/deltanufft;
  cout << " visi2dtacx/Info:   frequency bands ... " <<endl;
  cout << " Band 1: "<<freq1<<" +/-"<<deltafreq<<" MHz -> "<<JFmin<<" <=NumFreq<= "<<JFmax<<endl;
  cout << " Band 2: "<<freq2<<" +/-"<<deltafreq<<" MHz -> "<<JFminB<<" <=NumFreq<= "<<JFmaxB<<endl;

  // Bande de frequence de 2 MHz autour du 21 cm (1420 MHz -> JF=2792) 
  sa_size_t JFmin21=2776,JFmax21=2808;   // +/- 1 MHz autour de 1420 MHz 
  sa_size_t JFmin21_5=2710,JFmax21_5=2874;  //  +/- 5 MHz autour de 1420 MHz 
  cout << " pjHI +/-1 MHz @1420 MHz " << JFmin21<<" <=NumFreq<= "<<JFmax21
       << " pjHI +/-5 MHz @1420 MHz " << JFmin21_5<<" <=NumFreq<= "<<JFmax21_5 << endl;
 
  HiStatsInitiator _inia;
  //   TArrayInitiator  _inia;

  int rc = 0;
  try {
    ResourceUsage resu;

    //  Numero des lignes des auto-correlations dans la matrice des visibilites 
    sa_size_t KVAC[8]={0,8,15,21,26,30,33,35};
    //  Numero des lignes des cross-correlations 1H-2H 1H-3H 1H-4H 2H-3H 2H-4H 3H-4H dans la matrice des visibilites 
    sa_size_t KVCXHH[6]={1,2,3,9,10,16}; 

    cout << " visi2dtacx/Info: reading input gain file"<<gainfile<<" ..."; 
    Matrix gains;
    Vector gn;
    {
      PInPersist pin(gainfile);
      pin>>PPFNameTag("gains")>>gains;
      pin>>PPFNameTag("gn")>>gn;
      for (int p=0; p<8; p++) cout<<"gn["<<p<<"]="<<gn(p)<<" ";  
      // on transforme gains en facteur multiplicatif
      for(sa_size_t r=0;r<gains.NRows();r++) {
	for(sa_size_t c=0; c<gains.NCols(); c++) gains(r,c)=1./gains(r,c);
      }
      cout<<endl;
    }
    // Gains pour les cross-correlations g_ij(nu) = sqrt(g_ii(nu) 8 g_jj(nu) )
    Matrix gaincx(6, gains.NCols());
    for(sa_size_t c=0; c<gains.NCols(); c++)   {
      gaincx(0,c)=sqrt(gains(0,c)*gains(1,c));     //  1H-2H   
      gaincx(1,c)=sqrt(gains(0,c)*gains(2,c));     //  1H-3H   
      gaincx(2,c)=sqrt(gains(0,c)*gains(3,c));     //  1H-4H   
      gaincx(3,c)=sqrt(gains(1,c)*gains(2,c));     //  2H-3H   
      gaincx(4,c)=sqrt(gains(1,c)*gains(3,c));     //  2H-4H 
      gaincx(5,c)=sqrt(gains(2,c)*gains(3,c));     //  3H-4H     
    }
    cout << " DONE"<<endl;

    const char* dtnames[53]={"k","day","time",
			     "p1","p2","p3","p4","p5","p6","p7","p8",
			     "p1B","p2B","p3B","p4B","p5B","p6B","p7B","p8B",
			     "p1HI","p2HI","p3HI","p4HI","p5HI","p6HI","p7HI","p8HI",
			     "p1HI5","p2HI5","p3HI5","p4HI5","p5HI5","p6HI5","p7HI5","p8HI5",
                             "xh1h2","xh1h3","xh1h4","xh2h3","xh2h4","xh3h4",
                             "xh1h2B","xh1h3B","xh1h4B","xh2h3B","xh2h4B","xh3h4B",
			     "xh1h2HI","xh1h3HI","xh1h4HI","xh2h3HI","xh2h4HI","xh3h4HI"};

    DataTable dt(384);
    dt.AddIntegerColumn(dtnames[0]);
    dt.AddIntegerColumn(dtnames[1]);
    dt.AddFloatColumn(dtnames[2]);
    for(int kk=0; kk<32; kk++) 
      dt.AddFloatColumn(dtnames[3+kk]);
    for(int kk=0; kk<18; kk++) 
      dt.AddComplexColumn(dtnames[3+32+kk]);

    Range freqs(JFmin,JFmax);
    Range freqs21(JFmin21,JFmax21);
    Range freqs21_5(JFmin21_5,JFmax21_5);
    Range freqsB(JFminB,JFmaxB);

    vector<string> paths; 
    paths.push_back(path5); 
    paths.push_back(path6);
    VisiP4Reader vreader(paths, Imin,Imax,Istep,true);
    vreader.setPrintLevel(prtlev);
    bool fgok=true;
    TMatrix< complex<r_4> > vismtx;
    TMatrix< complex<r_4> > acsum;
    TMatrix< complex<r_4> > cxsum;
    double acdt[32];   // les 4*8=32 valeurs d'autocorrelation pour remplissage dans la table 
    complex<double> cxdt[18];   // les 3*6=18 valeurs de cross-correlations pour remplissage dans la table 

212
    TimeStamp dateobs, cfdate;
213 214 215 216
    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; 
217 218 219
    //----- 6 TimeFrequency maps 
    vector< TArray< complex<r_4> > > vtfm;
    
220 221
    // Getting en empty row_ptr
    DataTableRowPtr rowp = dt.EmptyRowPtr();
222 223 224 225 226

    //---- for the time-freqency map filling    
    sa_size_t TFMtmidx=0;
    sa_size_t tfmSX, tfmSY;

227
    while (fgok) {
228
      fgok=vreader.ReadNext(vismtx, cfdate, mttag);
229 230 231 232
      if (fgok)  {
	if (cnt==0)  {    //resizing matrices for sum of auto-correlations and sum of 6 cross-correlations 
	  acsum.SetSize(8, vismtx.NCols());
	  cxsum.SetSize(6, vismtx.NCols());
233 234 235 236 237 238
	  if (FgTFM) {  //allocating time-frequency maps 
	    tfmSX=(Imax-Imin)/Istep/deltaIavg;
	    tfmSY=vismtx.NCols()/TFMfbin;
	    cout << " visi2dtacx/Info: allocating Time-Freqncy 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) ); 
	  }
239 240
	}
	if (I==0) {   // start filling a new DataTable row 
241
	  dateobs=cfdate;
242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267
	  if (prtlev>0) 
	    cout << "visi2dtacx/Info:  dateobs="<<dateobs<<" SecondsPart()="<<dateobs.SecondsPart()<<endl;
	  for(int k=0; k<32; k++)  acdt[k]=0.;
	  for(int k=0; k<18; k++)  cxdt[k]=complex<r_8>(0.,0.);
	  acsum = complex<r_4>(0.,0.);
	  cxsum = complex<r_4>(0.,0.);
	}

	//   sum (integration) along the time axis 
	for(int k=0; k<8; k++)    acsum.Row(k) += vismtx.Row(KVAC[k]);     // Les auto-correlations 
	for(int k=0; k<6; k++)    cxsum.Row(k) += vismtx.Row(KVCXHH[k]);   // les cross-correlations 
	I++;    // incrementing DeltaTime counter 

	if (I==deltaIavg) {
	  //---- On s'occupe d'abord des autocorrelations P1 ... P8 
	  Vector vac(acsum.NCols());
	  for(int k=0; k<8; k++)  {
	  // correct for gain(nu)
	    for(sa_size_t c=0; c<vac.Size(); c++) vac(c)=acsum(k,c).real()*gains(k,c);
	    acdt[k] += vac(freqs).Sum();    // integration en bande 1
	    acdt[k+8] += vac(freqsB).Sum();    // integration en bande large B
	    acdt[k+16] += vac(freqs21).Sum();   // integration a +/- 1 MHz @ 1420 MHz 
	    acdt[k+24] += vac(freqs21_5).Sum();   // integration a +/- 5 MHz @ 1420 MHz 
	  }
	  //---- On s'occupe des cross-correlations  1H-2H ... 3H-4H 
	  TVector< complex<r_8> > vcx(acsum.NCols());
268
	  for(int k=0; k<6; k++)  {   // loop over the 6 Xcor 
269 270 271
	    for(sa_size_t c=0; c<vcx.Size(); c++) 
	      vcx(c)=complex<r_8>((r_8)cxsum(k,c).real(), (r_8)cxsum(k,c).imag())*complex<r_8>(gaincx(k,c), 0.);
	    cxdt[k] += vcx(freqs).Sum();    // integration en bande 1
272 273
	    cxdt[k+6] += vcx(freqsB).Sum();    // integration en bande large B
	    cxdt[k+12] += vcx(freqs21).Sum();   // integration a +/- 1 MHz @ 1420 MHz 
274 275 276 277 278 279 280 281 282 283 284

	    if (FgTFM && (TFMtmidx<tfmSX)) {  //filling time-frequency maps 
	      TArray< complex<r_4> > & tfmap = vtfm[k];
	      for(sa_size_t jy=0; jy<tfmSY; jy++) {
		complex<r_8> zz = vcx( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
		tfmap(TFMtmidx, jy) = complex<r_4>(zz.real(), zz.imag());
	      } 
	    }  //----- end of time-frequency maps filling 
	  }  //----- end of loop over the 6 Xcor 
	  TFMtmidx++;

285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300
	  // Moyennes dans la bande 1
	  double fnorm=(1./(double)deltaIavg)/(double)(JFmax-JFmin+1);
	  for(int k=0; k<8; k++)  acdt[k]*=fnorm; 
	  for(int k=0; k<6; k++)  cxdt[k]*=complex<r_8>(fnorm,0.);   
	  // Moyennes dans la bande 2 
	  fnorm=(1./(double)deltaIavg)/(double)(JFmaxB-JFminB+1);
	  for(int k=0; k<8; k++)  acdt[k+8]*=fnorm;
	  for(int k=0; k<6; k++)  cxdt[k+6]*=complex<r_8>(fnorm,0.);   
	  // Moyennes ds +/5 MHz @1420 , exclu +/- 1 MHz centre a 1420 
	  fnorm=(1./(double)deltaIavg)/(double)((JFmax21_5-JFmin21_5)-(JFmax21-JFmin21));
	  for(int k=0; k<8; k++)  acdt[k+24] = (acdt[k+24]-acdt[k+16])*fnorm;
	  // Moyennes ds +/1 MHz @1420 
	  fnorm=(1./(double)deltaIavg)/(double)(JFmax21-JFmin21+1);
	  for(int k=0; k<8; k++)  acdt[k+16]*=fnorm;
	  for(int k=0; k<6; k++)  cxdt[k+12]*=complex<r_8>(fnorm,0.);   
	  // remplissage datatable 
301 302 303 304
	  dt.NextRowPtr(rowp);
	  rowp(0)=(int_4)cnt;   
	  rowp(1)=(int_4)(dateobs.DaysPart()-dateorg.DaysPart());   
	  rowp(2)=(r_4)dateobs.SecondsPart();
305 306 307 308 309 310 311
	  for(int k=0; k<32; k++)   rowp(3+k) = (float)acdt[k];
	  for(int k=0; k<18; k++)   rowp(35+k) = complex<float>((float)cxdt[k].real(), (float)cxdt[k].imag());
	  //  ... done 
	  I=0;  cntnt++;
	}
	cnt++;
	if ((cntnt>0)&&(cntnt%10==0)&&(cntnt>pcntnt)) {
312 313
	  cout << " visi2dtacx/Info: fill-count="<<cntnt<<" fileCount="<<cnt<<" /Max="<<Imax
	       <<" DateObs="<<dateobs<<endl;
314 315 316 317 318 319 320
	  pcntnt=cntnt;
	}
      } 
    }
    cout << " visi2dtacx/Info: count="<<cnt<<" visimtx read "<<endl;
    dt.SetShowMinMaxFlag(true);
    cout << dt;
321
    cout<<"  visi2dtacx/Info: Saving auto-corr=f(time) datatable to PPF file "<<outfile<<endl;
322 323
    POutPersist po(outfile);
    po<<dt; 
324 325 326 327 328 329 330 331 332 333 334

    if (FgTFM) {   // --- renormalizing and saving time-frequency maps 
      cout<<"  visi2dtacx/Info: Saving time-frequency maps to PPF file "<<TFMoutfile<<endl;
      POutPersist potfm(TFMoutfile);
      const char* tfm_names[6]={"TFM_1H2H", "TFM_1H3H", "TFM_1H4H", "TFM_2H3H", "TFM_2H4H", "TFM_3H4H"};
      for(int k=0; k<6; k++)  {   // loop over the 6 Xcor 
	TArray< complex<r_4> > & tfmap = vtfm[k];
	tfmap /= complex<r_4>((r_4)(1./(double)deltaIavg/(double)TFMfbin), 0.);
	potfm << PPFNameTag(tfm_names[k]) << tfmap;
      }
    }
335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358
    //    resu.Update();
    cout << resu;   // Update est fait lors du print
  }
  catch (PException& exc) {
    cerr << " visi2dtacx.cc catched PException " << exc.Msg() << endl;
    rc = 77;
  }  
  catch (std::exception& sex) {
    cerr << "\n visi2dtacx.cc std::exception :" 
         << (string)typeid(sex).name() << "\n msg= " 
         << sex.what() << endl;
    rc = 78;
  }
  catch (...) {
    cerr << " visi2dtacx.cc catched unknown (...) exception  " << endl; 
    rc = 79; 
  } 

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

}