rdthermrfi.cc 6.61 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
// Utilisation de SOPHYA pour faciliter les tests ...
#include "sopnamsp.h"
#include "machdefs.h"

/* ---------------------------------------------------------- 
   Projet BAORadio/PAON4 - (C) LAL/IRFU  2008-2018

   Programme de lecture des fichiers matrices de visibilites de PAON4 
   produits par mfacq pour les voies Thermometres et RFI (9,10)
   R. Ansari  Univ. Paris Sud, & LAL CNRS/IN2P3 
   Juillet 2018 
   ---------------------------------------------------------- */

// include standard c/c++
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <iostream>
#include <string>

#include "pexceptions.h"
#include "array.h"
#include "histats.h"
#include "fitsioserver.h"

#include "resusage.h"

// include lecteur de fichiers visibilites 
#include "p4autils.h"
#include "visp4winreader.h"

int Usage(void)
{
  cout<<"--- rdthermrfi : Read PPF files for thermometer & RFI signal channnels and produce TF-Maps \n"<<endl;
  cout<<"  Usage: rdthermrfi [-options] \n";
  cout<<"  Specify -inth PathForThermometer or in datacards " <<endl;
  P4AnaParams::UsageOptions();
  cout<<endl;
  return 1;
}

//----------------------------------------------------
int main(int narg, const char* arg[])
{
  // --- Decoding parameters 
  if ((narg<2)||((narg>1)&&(strcmp(arg[1],"-h")==0)))  return Usage();
  P4AnaParams params;
  params.DecodeArgs(narg, arg);
  if (params.inpath_thermRFI_.length()<1) {
    cout<<"rdthermrfi/ERROR no path specified for thermometer/RFI signals ... (rdthermrfi -h for usage)" <<endl;
    return 2;
  }
  string outfile = params.outfile_;
  if (outfile.length()<1)  outfile = "thrfiavg.ppf";

  string fitsoutfile = params.fitsoutfile_;
  if (fitsoutfile.length()>=1) {
    fitsoutfile = "!"+fitsoutfile ; // adds '!' ?
  }

  int deltaIavg = params.TFMtimebin_;
  sa_size_t TFMfbin = params.TFMfreqbin_;
  
  int rc = 0;
  try {
    FitsIOServerInit();
    HiStatsInitiator _inia;

    ResourceUsage resu;

    VisiP4WindowReader wreader(params, true);  // true -> Read Therometer/RFI 

    long Imin = wreader.getReader().getSerialFirst();
    long Imax =  wreader.getReader().getSerialLast();
    long Istep =  wreader.getReader().getSerialStep();
    cout << "vis2ra/Info: processing visibility matrix serial/sequence number range "
	 <<Imin<<" <= seq <= " << Imax << " with step="<<Istep<<endl;
    cout << " WindowSize="<<wreader.getWindowSize()<<"  -> TotalNbWindows="<<wreader.getTotalNbWindows()<<endl;


    //vismtx: visib matrice, meanmtx: mean FFT coef, vismtx_center: visib matrix corrigee de l'offset (normalisee avec NPAQSUM)
    // un vecteur avec les temps 
    TVector< double > timevec(wreader.getTotalNbWindows()/deltaIavg); 
    TVector< double > ravec(wreader.getTotalNbWindows()/deltaIavg); 
    TMatrix< complex<r_4> > meanvismtx;
    TMatrix< complex<r_4> > visum; 

    TimeStamp dateobs, cfdate,datestart;
    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; 

    //----- 2 auto-corr (Thermometer/RFI) TimeFrequency maps
    TArray<r_4> TFM_therm;
    TArray<r_4> TFM_RFI;
    //----  1 cross-correlation TFM   (Thermometer <> RFI)
    TArray< complex<r_4> > TFM_cross;

    //---- for the time-freqency map filling    
    sa_size_t TFMtmidx=0;
    sa_size_t tfmSX, tfmSY;

    bool fgok=true;
    while (fgok) {
      //reads next visibility matrix window 
      fgok = wreader.Shift();
      if (!fgok)  break;
      
      TMatrix< complex<r_4> > vismtx = wreader.getAverageVisMtx(cfdate);
      if (cnt==0)  {    //resizing matrices for sum of auto-correlations and sum of 6 cross-correlations
	visum = vismtx;  meanvismtx=vismtx;
	tfmSX=wreader.getTotalNbWindows()/deltaIavg;
	tfmSY=vismtx.NCols()/TFMfbin;
	TFM_therm.SetSize2D(tfmSX, tfmSY);
	TFM_RFI.SetSize2D(tfmSX, tfmSY);
	TFM_cross.SetSize2D(tfmSX, tfmSY);
	dateobs=cfdate;
	continue;
      }
      if (I==0)  dateobs=cfdate; // start filling a new time bin
      // Accumulating 
      visum += vismtx;
      meanvismtx += vismtx;
      
      I++;    // incrementing DeltaTime counter 
      // we check that our time index did not go beyond the allocated array size (might not be necessary)
      if ((I==deltaIavg)&&(TFMtmidx>=tfmSX)) {  // Cela ne devrait pas arriver en principe 
	TFMtmidx++;	I=0;  
	cout << "rdthermrfi/Warning: something wrong in the logic , (TFMtmidx="<<TFMtmidx<<") >= (tfmSX="<<tfmSX<<")"
	     << " for read count="<<cnt<<endl; 
      }
      else if (I==deltaIavg) {   // Filling TFM maps
	TVector<r_4> vac0 = real(visum.Row(0));
	TVector<r_4> vac1 = real(visum.Row(1));
	TVector< complex<r_4> > vacX = visum.Row(2);
	for(sa_size_t jy=0; jy<tfmSY; jy++) {  // frequency binning 
	  TFM_therm(TFMtmidx, jy) = vac0( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
	  TFM_RFI(TFMtmidx, jy) = vac1( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
	  TFM_cross(TFMtmidx, jy) = vacX( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
	}
	double tdif =  cfdate.TimeDifferenceSeconds(cfdate,dateobs)/2.;
	timevec(TFMtmidx) = dateobs.TimeDifferenceSeconds(dateobs.ShiftSeconds (tdif ),datestart);	// centre du bin 
	ravec(TFMtmidx) = P4Coords::RAFromTimeTU(dateobs.ShiftSeconds (tdif ));
	TFMtmidx++;
	//  ... done 
	I=0;  cntnt++;
	visum = complex<r_4>(0.,0.);
      }
      cnt++;
      if ((cnt>0)&&(cntnt%10==0)&&(cntnt>pcntnt)) {
	cout<<"rdthermrfi/Info: TFM-Map fill cnt="<<cntnt<<" VisMtxCount="<<cnt
	    <<" /Max="<<wreader.getTotalNbWindows()<<" DateObs="<<dateobs<<endl;
	pcntnt=cntnt;
      }
    }  // End of redaing visibility matrices loop
    
    cout<<"rdthermrfi/Info: count="<<cnt*wreader.getWindowSize()<<" Visibility Matrices read "<<endl;
    if (cnt>0)  { 
      meanvismtx /= complex<r_4>((r_4)cnt, (r_4)0.);
      cout<<"  rdthermrfi/Info: Saving vismtx_mean to PPF file "<<outfile<<endl;
      POutPersist po(outfile);
      po<<PPFNameTag("meanvis_Th_RFI")<<meanvismtx;
      cout<<"  rdthermrfi/Info: TFM_therm , TFM_RFI, TFM_cross (TFM_ThXRFI) to PPF file "<<outfile<<endl;
      po<<PPFNameTag("TFM_therm")<<TFM_therm;
      po<<PPFNameTag("TFM_RFI")<<TFM_RFI;
      po<<PPFNameTag("TFM_ThXRFI")<<TFM_cross;
    } 
    //    resu.Update();
    cout << resu;   // Update est fait lors du print
  }
  catch (PException& exc) {
    cerr << " rdthermrfi.cc catched PException " << exc.Msg() << endl;
    rc = 77;
  }  
  catch (std::exception& sex) {
    cerr << "\n rdthermrfi.cc std::exception :" 
         << (string)typeid(sex).name() << "\n msg= " 
         << sex.what() << endl;
    rc = 78;
  }
  catch (...) {
    cerr << " rdthermrfi.cc catched unknown (...) exception  " << endl; 
    rc = 79; 
  } 

  cout << ">>>> rdthermrfi.cc ------- END ----------- RC=" << rc << endl;
  return rc;

}