visi2ntac.cc 6.43 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 

   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"

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

/* --Fonction-- */
int Usage(bool fgea)
{
  cout << " --- visi2ntac.cc : Read PPF files produced by mfacq time-frequency\n" << endl;
  if (fgea) cout << " visi2ntac: Missing or wrong argument ! " << endl;
51
  cout << " Usage: visi2ntac InPathBAO5 InPathBAO6 Imin,Imax OutPPFFile [Jf1,Jf2] [DeltaIAvg] [PrtLev=0] \n"
52 53 54
       << "  InPathBAO5/6: path for input vismtx_J_iii.ppf on bao5/6 \n"
       << "  Imin,Imax : read files for Imin<=iii<=Imax  iii+=step \n"
       << "  OutPPFFile: Output PPF file name \n"
55
       << "  Jf1,Jf2: frequency index (Jf1<=f<=Jf2) range for power integration, def=768,2304 \n"
56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78
       << "  DeltaIAvg: compute average power every DeltaI input vismtx files, def=1 \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<5) 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 outfile=arg[4];
  // frequency range definition 
79
  sa_size_t JFmin=768,JFmax=2304;
80
  int jf1,jf2;
81
  if ((narg>5)&&(strcmp(arg[5],".")!=0))  sscanf(arg[5],"%d,%d",&jf1,&jf2);
82 83 84 85 86
  JFmin=jf1;  JFmax=jf2;
  // time averaging interval definition, in term of visibility files
  int deltaIavg=1;
  if (narg>6)  deltaIavg=atoi(arg[6]);
  if (deltaIavg<1)  deltaIavg=1;
87 88 89
  int prtlev=1;
  if (narg>7) prtlev=atoi(arg[7]);

90 91
  cout << " visi2ntac/Info: Path BAO5:"<<path5<<" BAO6:"<<path6<<"\n"
       << " Imin,max,step="<<Imin<<","<<Imax<<","<<Istep<<"  OutFile:"<<outfile<<"\n"
92
       << JFmin<<" <=NumFreq<= "<<JFmax<<"  DeltaIAvg="<<deltaIavg<<" PrtLev="<<prtlev<<endl;
93 94 95 96 97 98 99

  // Bande de frequence de 2 MHz autour du 21 cm (1420 MHz -> JF=2758) 
  sa_size_t JFmin21=2748,JFmax21=2768;   // +/- 1 MHz autour de 1420 MHz 
  sa_size_t JFmin21_5=2676,JFmax21_5=2840;  //  +/- 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;
 
100 101 102 103 104 105 106 107 108 109
  HiStatsInitiator _inia;
  //   TArrayInitiator  _inia;

  int rc = 0;
  try {
    ResourceUsage resu;
    vector<string> paths; 
    paths.push_back(path5); 
    paths.push_back(path6);

110 111 112 113 114 115 116
    const char* ntnames[26]={"k","time",
			     "p1","p2","p3","p4","p5","p6","p7","p8",
			     "p1HI","p2HI","p3HI","p4HI","p5HI","p6HI","p7HI","p8HI",
			     "p1HI5","p2HI5","p3HI5","p4HI5","p5HI5","p6HI5","p7HI5","p8HI5"};

    NTuple nt(26,ntnames);
    r_8 xnt[30];
117 118 119
    //  Numero de ligne des auto-correlations dans la matrice des visibilites 
    sa_size_t KVAC[8]={0,8,15,21,26,30,33,35};
    Range freqs(JFmin,JFmax);
120 121 122
    Range freqs21(JFmin21,JFmax21);
    Range freqs21_5(JFmin21_5,JFmax21_5);

123
    VisiP4Reader vreader(paths, Imin,Imax,Istep,true);
124
    vreader.setPrintLevel(prtlev);
125 126 127 128
    bool fgok=true;
    TMatrix< complex<r_4> > vismtx;
    TimeStamp dateobs;
    double mttag;
129
    int cnt=0, cntnt=0, pcntnt=0;
130 131 132 133 134 135
    int I=0; 
    while (fgok) {
      fgok=vreader.ReadNext(vismtx, dateobs, mttag);
      if (fgok)  {
	if (I==0) {
	  xnt[0]=cnt;   xnt[1]=dateobs.SecondsPart();
136 137
	  if (prtlev>0) 
	    cout << "visi2ntac/Info:  dateobs="<<dateobs<<" SecondsPart()="<<dateobs.SecondsPart()<<endl;
138
	  for(int k=0; k<24; k++)  xnt[2+k]=0.;
139
	}
140
	for(int k=0; k<8; k++)  {     // integration en bande large 
141 142 143
	  TVector< complex<r_4> > vac = vismtx.Row(KVAC[k]);
	  xnt[2+k] += vac(freqs).Sum().real();
	}
144 145 146 147 148 149 150 151
	for(int k=0; k<8; k++)  {     // integration a +/- 1 MHz @ 1420 MHz 
	  TVector< complex<r_4> > vac = vismtx.Row(KVAC[k]);
	  xnt[2+8+k] += vac(freqs21).Sum().real();
	}
	for(int k=0; k<8; k++)  {     // integration a +/- 5 MHz @ 1420 MHz 
	  TVector< complex<r_4> > vac = vismtx.Row(KVAC[k]);
	  xnt[2+16+k] += vac(freqs21_5).Sum().real();
	}
152 153 154 155
	I++; 
	if (I==deltaIavg) {
	  double fnorm=(1./(double)deltaIavg)/(double)(JFmax-JFmin+1);
	  for(int k=0; k<8; k++)  xnt[2+k]*=fnorm;
156 157 158 159
	  fnorm=(1./(double)deltaIavg)/(double)(JFmax21-JFmin21+1);
	  for(int k=0; k<8; k++)  xnt[2+8+k]*=fnorm;
	  fnorm=(1./(double)deltaIavg)/(double)(JFmax21_5-JFmin21_5+1);
	  for(int k=0; k<8; k++)  xnt[2+16+k]*=fnorm;
160
	  nt.Fill(xnt);
161
	  I=0;  cntnt++;
162 163
	}
	cnt++;
164 165 166 167
	if ((cntnt>0)&&(cntnt%10==0)&&(cntnt>pcntnt)) {
	  cout << " visi2ntac/Info: fill-count="<<cntnt<<" DateObs="<<dateobs<<endl;
	  pcntnt=cntnt;
	}
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
      } 
    }
    cout << " visi2ntac/Info: count="<<cnt<<" visimtx read "<<endl;
    cout << nt;
    cout<<"  rdvisip4/Info: Saving auto-corr=f(time) ntuple to PPF file "<<outfile<<endl;
    POutPersist po(outfile);
    po<<nt; 
    //    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;

}