rdthermrfi.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 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
// 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();
78
    cout << "rdthermrfi/Info: processing visibility matrix serial/sequence number range "
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
	 <<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);
113 114

      if (cnt==0)  {    //resizing matrices for sum of auto-correlations and sum of  cross-correlations
115 116 117 118 119 120 121
	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;
122 123
	// recupere le jour de depart @ 0h
	datestart = TimeStamp(cfdate.DaysPart(),0.);
124 125 126 127 128
	// y avait cnt++ & continue ?
      }// end of resizing 
      else{
	visum += vismtx;
	meanvismtx += vismtx;
129 130 131
      }
      if (I==0)  dateobs=cfdate; // start filling a new time bin
      // Accumulating 
132
      
133 134 135 136 137 138 139 140 141
      
      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
142
	TVector<r_4> vac0 = real(visum.Row(0));   // Voie 1 : Thermometre 
perdereau's avatar
perdereau committed
143 144
	TVector<r_4> vac1 = real(visum.Row(2));   // Voie 2 : RFI
	TVector< complex<r_4> > vacX = visum.Row(1);  //  Cx Thermometere X RFI 
145 146 147 148 149 150 151 152 153 154 155
	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++;
156
	visum = complex<r_4>(0.,0.); // RAZ
157 158 159 160 161 162 163
      }
      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;
      }
164
    }  // End of reading visibility matrices loop
165 166 167 168 169 170 171 172 173 174 175
    
    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;
perdereau's avatar
perdereau committed
176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213
      po << PPFNameTag("TimeVec") << timevec ;
      po << PPFNameTag("RAVec") << ravec ;

      P4FreqBand myp4fre;
      TVector <double> avg_freqs( myp4fre.getP4NbFreqChannels()/TFMfbin);
      double frbase =  myp4fre.freqstart_ + myp4fre.getP4FreqResolution()/2. ;
      for (int kf=0 ; kf< myp4fre.getP4NbFreqChannels()/TFMfbin ; kf++,frbase += myp4fre.getP4FreqResolution()*TFMfbin )
	avg_freqs(kf) = frbase ;
    
      po << PPFNameTag("FreqVec") << avg_freqs;

      FitsInOutFile  * fos = NULL ;
      
      if (fitsoutfile.length()>=1){
	vector <string> ext_names;

	cout << " fitsoutfile :" <<fitsoutfile<<":"<< endl;
	fos = new FitsInOutFile(fitsoutfile, FitsInOutFile::Fits_Create);

	ext_names.push_back("TFM_therm");
	(*fos) << TFM_therm;
	ext_names.push_back("TFM_RFI");
	(*fos) << TFM_RFI;
	ext_names.push_back("TFM_ThXRFI_R");
	(*fos) << real(TFM_cross);
	ext_names.push_back("TFM_ThXRFI_I");
	(*fos) << imag(TFM_cross);

	ext_names.push_back("Frequences");
	(*fos)<< avg_freqs ;
	ext_names.push_back("RAs");
	(*fos)<<  ravec;
	ext_names.push_back("Times");
	(*fos)<<  timevec;
	cout << " number of objs in fits "<< ext_names.size() << endl;
	cout << ext_names << endl;
	delete(fos);
      }
214
    }// end if cnt>0 
215 216
    //    resu.Update();
    cout << resu;   // Update est fait lors du print
217
  }//end try
218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237
  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;

}