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

rdthermrfi.cc 7.88 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
// 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;
121
	cnt++;
122 123 124 125 126 127 128 129 130 131 132 133 134 135 136
	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
137
	TVector<r_4> vac0 = real(visum.Row(0));   // Voie 1 : Thermometre 
perdereau's avatar
perdereau committed
138 139
	TVector<r_4> vac1 = real(visum.Row(2));   // Voie 2 : RFI
	TVector< complex<r_4> > vacX = visum.Row(1);  //  Cx Thermometere X RFI 
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
	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;
perdereau's avatar
perdereau committed
171 172 173 174 175 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
      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);

      }


212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235
    } 
    //    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;

}