Docker-in-Docker (DinD) capabilities of public runners deactivated. More info

visi2ntac.cc 7.46 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 GainPPFFile OutPPFFile Jf1A,Jf2A JF1B,Jf2B DeltaIAvg [PrtLev=0] \n"
52 53
       << "  InPathBAO5/6: path for input vismtx_J_iii.ppf on bao5/6 \n"
       << "  Imin,Imax : read files for Imin<=iii<=Imax  iii+=step \n"
54
       << "  GainPPFFile: Input gains PPF file name \n"
55
       << "  OutPPFFile: Output PPF file name \n"
56
       << "  Jf1A/B,Jf2A/B : frequency index (Jf1<=f<=Jf2) range for power integration, A:1300,1460 B:2130,2290 \n"
57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
       << "  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);
72
  if (narg<9) return Usage(true);
73 74 75 76 77
  string path5 = arg[1];
  string path6 = arg[2];
  int Imin,Imax,Istep=1; 
  sscanf(arg[3],"%d,%d",&Imin,&Imax);
  Istep=1;
78 79
  string gainfile=arg[4];
  string outfile=arg[5];
80
  // frequency range definition 
81
  sa_size_t JFmin=1300,JFmax=1460;
82
  int jf1,jf2;
83
  if ((narg>6)&&(strcmp(arg[6],"-")!=0))  sscanf(arg[6],"%d,%d",&jf1,&jf2);
84
  JFmin=jf1;  JFmax=jf2;
85 86 87
  sa_size_t JFminB=2130,JFmaxB=2290;  
  if ((narg>7)&&(strcmp(arg[7],"-")!=0))  sscanf(arg[7],"%d,%d",&jf1,&jf2);
  JFminB=jf1;  JFmaxB=jf2;
88
  // time averaging interval definition, in term of visibility files
89 90
  int deltaIavg=10;
  if (narg>8)  deltaIavg=atoi(arg[8]);
91
  if (deltaIavg<1)  deltaIavg=1;
92
  int prtlev=1;
93
  if (narg>9) prtlev=atoi(arg[9]);
94

95 96
  cout << " visi2ntac/Info: Path BAO5:"<<path5<<" BAO6:"<<path6<<"\n"
       << " Imin,max,step="<<Imin<<","<<Imax<<","<<Istep<<"  OutFile:"<<outfile<<"\n"
97 98 99 100 101 102 103 104 105 106 107 108 109
       << JFmin<<" <=NumFreq<= "<<JFmax<<"  "<<JFminB<<" <=NumFreqB<= "<<JFmaxB
       <<"  DeltaIAvg="<<deltaIavg<<" PrtLev="<<prtlev<<endl;

  cout << " visi2ntac/Info: reading input gain file"<<gainfile<<endl;
  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)<<" ";  
    cout<<endl;
  }
110 111 112
  // 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 
113 114 115
  cout << " pjHI +/-1 MHz @1420 MHz " << JFmin21<<" <=NumFreq<= "<<JFmax21
       << " pjHI +/-5 MHz @1420 MHz " << JFmin21_5<<" <=NumFreq<= "<<JFmax21_5 << endl;
 
116 117 118 119 120 121 122 123 124 125
  HiStatsInitiator _inia;
  //   TArrayInitiator  _inia;

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

126
    const char* ntnames[34]={"k","time",
127 128
			     "p1","p2","p3","p4","p5","p6","p7","p8",
			     "p1HI","p2HI","p3HI","p4HI","p5HI","p6HI","p7HI","p8HI",
129 130
			     "p1HI5","p2HI5","p3HI5","p4HI5","p5HI5","p6HI5","p7HI5","p8HI5",
			     "p1B","p2B","p3B","p4B","p5B","p6B","p7B","p8B" };
131

132 133
    NTuple nt(34,ntnames);
    r_8 xnt[40];
134 135 136
    //  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);
137 138
    Range freqs21(JFmin21,JFmax21);
    Range freqs21_5(JFmin21_5,JFmax21_5);
139
    Range freqsB(JFminB,JFmaxB);
140

141
    VisiP4Reader vreader(paths, Imin,Imax,Istep,true);
142
    vreader.setPrintLevel(prtlev);
143 144 145 146
    bool fgok=true;
    TMatrix< complex<r_4> > vismtx;
    TimeStamp dateobs;
    double mttag;
147
    int cnt=0, cntnt=0, pcntnt=0;
148 149 150 151 152 153
    int I=0; 
    while (fgok) {
      fgok=vreader.ReadNext(vismtx, dateobs, mttag);
      if (fgok)  {
	if (I==0) {
	  xnt[0]=cnt;   xnt[1]=dateobs.SecondsPart();
154 155
	  if (prtlev>0) 
	    cout << "visi2ntac/Info:  dateobs="<<dateobs<<" SecondsPart()="<<dateobs.SecondsPart()<<endl;
156
	  for(int k=0; k<24; k++)  xnt[2+k]=0.;
157
	}
158
	for(int k=0; k<8; k++)  {    
159 160 161 162 163 164 165 166
	  TVector< complex<r_4> > vacz = vismtx.Row(KVAC[k]);
	  Vector vac(vismtx.NCols());
	  // correct for gain(nu)
	  for(sa_size_t c=0; c<vac.Size(); c++) vac(c)=vacz(c).real()/gains(k,c);
	  xnt[2+k] += vac(freqs).Sum();    // integration en bande large 
	  xnt[2+8+k] += vac(freqs21).Sum();   // integration a +/- 1 MHz @ 1420 MHz 
	  xnt[2+16+k] += vac(freqs21_5).Sum();   // integration a +/- 5 MHz @ 1420 MHz 
	  xnt[2+24+k] += vac(freqsB).Sum();    // integration en bande large B
167
	}
168 169
	I++; 
	if (I==deltaIavg) {
170
	  // puissance moyenne dans la bande large
171 172
	  double fnorm=(1./(double)deltaIavg)/(double)(JFmax-JFmin+1);
	  for(int k=0; k<8; k++)  xnt[2+k]*=fnorm;
173 174 175 176
	  // puissance moyenne 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++)  xnt[2+16+k] = (xnt[2+16+k]-xnt[2+8+k])*fnorm;
	  // puissance moyenne ds +/1 MHz @1420 
177 178
	  fnorm=(1./(double)deltaIavg)/(double)(JFmax21-JFmin21+1);
	  for(int k=0; k<8; k++)  xnt[2+8+k]*=fnorm;
179 180 181
	  // puissance moyenne en bande large B
	  fnorm=(1./(double)deltaIavg)/(double)(JFmaxB-JFminB+1);
	  for(int k=0; k<8; k++)  xnt[2+24+k]*=fnorm;
182
	  nt.Fill(xnt);
183
	  I=0;  cntnt++;
184 185
	}
	cnt++;
186 187 188 189
	if ((cntnt>0)&&(cntnt%10==0)&&(cntnt>pcntnt)) {
	  cout << " visi2ntac/Info: fill-count="<<cntnt<<" DateObs="<<dateobs<<endl;
	  pcntnt=cntnt;
	}
190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220
      } 
    }
    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;

}