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

visi2ntac.cc 8.05 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
// 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"
41 42
#include "p4autils.h"
#include "p4gnugain.h"
43 44 45 46 47 48 49 50 51 52

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

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

108 109
  P4gnuGain p4g(gainfile);
  Matrix const & gains = p4g.getGainMatrix();
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[35]={"k","day","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
    NTuple nt(35,ntnames,384,false);   // float numbers have enough precision for us
133
    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 142 143
    VisiP4ReaderBase * reader = VisiP4ReaderBase::getReader(paths);
    VisiP4ReaderBase & vreader = (*reader);
    vreader.setFreqReordering(fgreorderfreq);
144
    cout << " vreader.SelectSerialNum(Imin="<<Imin<<" ,Imax="<<Imax<<" ,Istep="<<Istep<<")"<<endl;
145
    vreader.SelectSerialNum(Imin,Imax,Istep);
146
    vreader.setPrintLevel(prtlev);
147
   
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;

}