visi2dtacx.cc 11.8 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 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 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 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 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
// 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;
  cout << " Usage: visi2dtacx InPathBAO5 InPathBAO6 Imin,Imax GainPPFFile OutPPFFile Freq1 Freq2 DeltaFreq DeltaIAvg [PrtLev=0] \n"
       << "  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);
  string path5 = arg[1];
  string path6 = arg[2];
  int Imin,Imax,Istep=1; 
  sscanf(arg[3],"%d,%d",&Imin,&Imax);
  Istep=1;
  string gainfile=arg[4];
  string outfile=arg[5];
  // frequency range definition 
  double deltafreq=1.;
  double freq1=1375.;
  double freq2=1410.;
  if ((narg>7)&&(strcmp(arg[6],"-")!=0))  freq1=atof(arg[6]);
  if ((narg>8)&&(strcmp(arg[7],"-")!=0))  freq2=atof(arg[7]);
  if ((narg>6)&&(strcmp(arg[8],"-")!=0))  deltafreq=atof(arg[8]);
  // time averaging interval definition, in term of visibility files
  int deltaIavg=10;
  if ((narg>9)&&(strcmp(arg[9],"-")!=0))  deltaIavg=atoi(arg[9]);
  if (deltaIavg<1)  deltaIavg=1;
  int prtlev=1;
  if (narg>10) prtlev=atoi(arg[10]);

  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;


  //---- 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 

195
    TimeStamp dateobs, cfdate;
196 197 198 199 200 201 202 203 204
    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; 

    // Getting en empty row_ptr
    DataTableRowPtr rowp = dt.EmptyRowPtr();
    // 
    while (fgok) {
205
      fgok=vreader.ReadNext(vismtx, cfdate, mttag);
206 207 208 209 210 211
      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());
	}
	if (I==0) {   // start filling a new DataTable row 
212
	  dateobs=cfdate;
213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242
	  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());
	  for(int k=0; k<6; k++)  {
	    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
243 244
	    cxdt[k+6] += vcx(freqsB).Sum();    // integration en bande large B
	    cxdt[k+12] += vcx(freqs21).Sum();   // integration a +/- 1 MHz @ 1420 MHz 
245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261
	  }
	  // 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 
262 263 264 265
	  dt.NextRowPtr(rowp);
	  rowp(0)=(int_4)cnt;   
	  rowp(1)=(int_4)(dateobs.DaysPart()-dateorg.DaysPart());   
	  rowp(2)=(r_4)dateobs.SecondsPart();
266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307
	  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)) {
	  cout << " visi2dtacx/Info: fill-count="<<cntnt<<" DateObs="<<dateobs<<endl;
	  pcntnt=cntnt;
	}
      } 
    }
    cout << " visi2dtacx/Info: count="<<cnt<<" visimtx read "<<endl;
    dt.SetShowMinMaxFlag(true);
    cout << dt;
    cout<<"  visi2dtacx/Info: Saving auto-corr=f(time) ntuple to PPF file "<<outfile<<endl;
    POutPersist po(outfile);
    po<<dt; 
    //    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;

}