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

visi2ntac.cc 8.04 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
       << JFmin<<" <=NumFreq<= "<<JFmax<<"  "<<JFminB<<" <=NumFreqB<= "<<JFmaxB
       <<"  DeltaIAvg="<<deltaIavg<<" PrtLev="<<prtlev<<endl;

  cout << " visi2ntac/Info: reading input gain file"<<gainfile<<endl;
  Matrix gains;
  Vector gn;
103
  double GAINTHR = 0.005;  // limite inferieur du gain avant faire 1/G
104 105 106 107 108
  {
    PInPersist pin(gainfile);
    pin>>PPFNameTag("gains")>>gains;
    pin>>PPFNameTag("gn")>>gn;
    for (int p=0; p<8; p++) cout<<"gn["<<p<<"]="<<gn(p)<<" ";  
109 110
    // on transforme gains en facteur multiplicatif
    for(sa_size_t r=0;r<gains.NRows();r++) {
111
      for(sa_size_t c=0; c<gains.NCols(); c++) gains(r,c)=((gains(r,c)>GAINTHR)?1./gains(r,c):0.);
112
    }
113 114
    cout<<endl;
  }
115 116 117
  // 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 
118 119 120
  cout << " pjHI +/-1 MHz @1420 MHz " << JFmin21<<" <=NumFreq<= "<<JFmax21
       << " pjHI +/-5 MHz @1420 MHz " << JFmin21_5<<" <=NumFreq<= "<<JFmax21_5 << endl;
 
121 122 123 124 125 126 127 128 129 130
  HiStatsInitiator _inia;
  //   TArrayInitiator  _inia;

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

131
    const char* ntnames[35]={"k","day","time",
132 133
			     "p1","p2","p3","p4","p5","p6","p7","p8",
			     "p1HI","p2HI","p3HI","p4HI","p5HI","p6HI","p7HI","p8HI",
134 135
			     "p1HI5","p2HI5","p3HI5","p4HI5","p5HI5","p6HI5","p7HI5","p8HI5",
			     "p1B","p2B","p3B","p4B","p5B","p6B","p7B","p8B" };
136

137
    NTuple nt(35,ntnames,384,false);   // float numbers have enough precision for us
138
    r_8 xnt[40];
139 140 141
    //  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);
142 143
    Range freqs21(JFmin21,JFmax21);
    Range freqs21_5(JFmin21_5,JFmax21_5);
144
    Range freqsB(JFminB,JFmaxB);
145

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

}