source: Sophya/trunk/AddOn/TAcq/brprocGain.cc

//----------------------------------------------------------------
// Projet BAORadio - (C) LAL/IRFU  2008-2010
// Classes de threads de traitement pour BAORadio 
//----------------------------------------------------------------

#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include <fstream>
#include <signal.h>
#include <algorithm>

#include "pexceptions.h"
#include "tvector.h"
#include "ntuple.h"
#include "datatable.h"
#include "histos.h"
#include "fioarr.h"
#include "matharr.h"
#include "timestamp.h"
#include "ctimer.h"
#include "fftpserver.h"
#include "fitsarrhand.h"

#include "FFTW/fftw3.h"


#include "pciewrap.h"
#include "brpaqu.h"

#include "brprocGain.h"




//---------------------------------------------------------------------
// Classe de traitement de spectres - 
// Calcul de spectres de gain
//---------------------------------------------------------------------
/*!
   \class BRGainCalculator
   \ingroup TAcq

   \brief Gain/normalization computation on spectra (developped for Nancay data)
*/
/* --Methode-- */
BRGainCalculator::BRGainCalculator(RAcqMemZoneMgr& memgr, string outpath, uint_4 nmean, 
                                           bool fgdatafft, bool fgsinglechan)
  : BRMeanSpecCalculator(memgr,outpath,nmean,fgdatafft,fgsinglechan), forceMedianFreqFilter_(false), medianFilterHalfWidth_(50), nbwin4mean_(nmean)
{

  //  cout << "(JEC): BRGainCalculator Ctor " << endl;

  
  setNameId("gainCalc",1);

  uint_4 nb_octets_entrop = 0; //this value is valid for Dec. 2010 data at Nancay
  
  const char* venvp = NULL; 
  venvp=getenv("BRANA_NBYTECUT");
  if (venvp!=NULL){
    nb_octets_entrop = atoi(venvp);
    cout << "BRGainCalculator : BRANA_NBYTECUT : " << nb_octets_entrop << endl;
  }

  BRPaquet paq(memgr_.PaqSize()-nb_octets_entrop);

  if (fgsinglechannel_) {
    medfiltspecmtx_.SetSize(memgr_.NbFibres(), paq.DataSize()/2);     
  }
  else {
    medfiltspecmtx_.SetSize(2*memgr_.NbFibres(), paq.DataSize()/4);     
  }
  //init
  medfiltspecmtx_ = (r_4)(0.);
  nummedianfile_ = 0;
  nbtot_specwin_ = 0;
  windowTime_ = TimeStamp(); //default

}


/* --Methode-- */
BRGainCalculator::~BRGainCalculator()
{
  cout << "(JEC): BRGainCalculator Dtor:  nbtot_specwin = " << nbtot_specwin_
       << "nbwin4mean = " << nbwin4mean_ << endl;
  if(nbtot_specwin_>nbwin4mean_*0.5)SaveMedianSpectra();
}



/* --Methode-- */
void BRGainCalculator::ProcSpecWin(uint_8 numpaqstart, uint_8 numpaqend)
{

  //   cout << "(JEC): BRGainCalculator ProcSpecWin: "<<   nbtot_specwin_ 
  //        << endl;
  


  //JEC 13/12/10 save the filtered spectra
  if ( (nbtot_specwin_>0)&&(nbtot_specwin_%nbwin4mean_==0) ) SaveMedianSpectra();

  // Get time of window
  windowTime_ = getObsTime();

 if (prtlev_>1)  {
      cout << " BRGainCalculator::ProcSpecWin() num_win=" << nbtot_specwin_ << " numpaqstart=" << numpaqstart 
           << " numpaqend=" << numpaqend  << " nbwin4mean = " << nbwin4mean_ << endl;
      cout << " (1) ObsTime=" << windowTime_ << " TimeTag=" << getCurTimeTagSeconds() << " s. FrameCounter=" 
           << getCurFrameCounter() << endl;
  }  
  //DBG  cout << "BRGainCalculator::ProcSpecWin()/Debug: numpaqstart=" << numpaqstart 
  //DBG     << " numpaqend=" << numpaqend << endl;


  //Implement the median filter at a fixed freq. but on all paquets of the window
  //then apply a median filter on the frequencies
  for(sa_size_t i=0; i<spec_window_.SizeY(); i++)  { //loop on channel (voie)
    TVector<r_4> medianFilterPaq(spec_window_.SizeX());
    //we keep freq. = 0 for compatibility
    for (sa_size_t f=0; f<spec_window_.SizeX(); f++) {
      std::vector<r_4> val1;
      TVector<r_4> tval1(spec_window_(Range(f),Range(i),Range::all()).CompactAllDimensions());
      tval1.FillTo(val1);
      std::nth_element(val1.begin(),
                       val1.begin()+val1.size()/2,
                       val1.end());
      medianFilterPaq(f) = *(val1.begin()+val1.size()/2);
    }//eo loop on frequencies



    if(forceMedianFreqFilter_) {
      //perform a median filter on the frequencies
      TVector<r_4> medianFilterFreq = medfiltspecmtx_.Row(i);
      sa_size_t hww = medianFilterHalfWidth_; //half freq window for filtering
      sa_size_t fMin = 0.;
      sa_size_t fMax = spec_window_.SizeX()-1;
      for (sa_size_t f=0; f<spec_window_.SizeX(); f++) {
        sa_size_t first = (f-hww >= fMin) ? f-hww : fMin;
        sa_size_t last = (f+hww <= fMax) ? f+hww : fMax;
        std::vector<r_4> val2;
        medianFilterPaq(Range(first,last)).FillTo(val2);
        std::nth_element(val2.begin(),
                         val2.begin()+val2.size()/2,
                         val2.end());
        medianFilterFreq(f) = *(val2.begin()+val2.size()/2); //Notice the "+" for later take the mean
      }//eo loop on frequencies
    } else {
      //      cout<< "(JEC) ProcSpecWin debug medfiltspecmtx/medianFilterPaq "<< endl;
      //      medfiltspecmtx_.Row(i).Show();
      //      medianFilterPaq.Transpose().Show();
      medfiltspecmtx_.Row(i) += medianFilterPaq.Transpose();
    }
  }//eo loop on channels

  //We always end here with a medfiltspecmtx Matrix Row=[0,NCha-1] x Col=[0,fMax]

  
  if (nbtot_specwin_<nmaxfiles_specw_)  SaveSpectraWindow();
  

  nbtot_specwin_++;
  return;
}


/* --Methode-- */
//JEC 13/12/10 Start
void BRGainCalculator::SaveMedianSpectra()
{

  cout << "(JEC): BRGainCalculator   SaveMedianSpectra ("<<nummedianfile_<<") nbwin = "<< nbwin4mean_  
       << " nbtot_specwin = " <<  nbtot_specwin_
       << endl;
  if (prtlev_>1)  {
    cout << " BRGainCalculator:SaveMedianSpectra() nbwin4mean = " << nbwin4mean_ << endl
         << " (2) ObsTime=" << windowTime_ << " Total Intensity " << 0.5*medfiltspecmtx_.Sum() << endl;
  }  
  
  //JEC compute power spectrum
  



  for(sa_size_t ir=0; ir< medfiltspecmtx_.NRows(); ir++){
    char buff[32];
    sprintf(buff,"NPAQSUM_%d",(int)ir);
    medfiltspecmtx_.Info()["NPAQSUM"] = nbtot_specwin_;
    medfiltspecmtx_.Info()[buff] = nbtot_specwin_;
    if ( nbtot_specwin_ > 0) {
      medfiltspecmtx_.Row(ir)  /= (r_4)nbtot_specwin_;
    }
  }//eo loop on channels
  //JEC timeStamp in the day @ millisec: faudra gerer le passage a minuit!
  //Reza : getObsTimeSeconds() fournit le temps ecoule en seconde a partir de 0h00 de la date du premier fichier
  medfiltspecmtx_.Info()["TIMEWIN"] = (windowTime_.SecondsPart())*1000.;
  if (nbpmoyttfc_>1)  {
    moyfc_/=(double)nbpmoyttfc_;
    moytt_/=(double)nbpmoyttfc_;
  }
  string ikey,ikdesc;    
  ikey="MeanFC";  ikdesc="Mean FrameCounter";
  medfiltspecmtx_.Info().SetI(ikey,moyfc_); 
  medfiltspecmtx_.Info().SetComment(ikey,ikdesc);
  ikey="MeanTT";  ikdesc="Mean TimeTag";
  medfiltspecmtx_.Info().SetD(ikey,moytt_);
  medfiltspecmtx_.Info().SetComment(ikey,ikdesc);

  char nfile[64];
  string flnm;
  {
    sprintf(nfile,"medfiltmtx%d.fits",nummedianfile_);
    flnm="!"+outpath_+nfile;
    FitsInOutFile fos(flnm,FitsInOutFile::Fits_Create);
    fos << medfiltspecmtx_;
    cout << nummedianfile_ << "-BRGainCalculator::ProcSpecWin() save filtered spectra matrix in " 
         << flnm << endl;
  }
  
  //increment the file index
  nummedianfile_++;

  if (dtpms_!=NULL)  FillPwrTmDTable(medfiltspecmtx_);
    
  //reset the matirx
  medfiltspecmtx_ = (r_4)(0.);
  //reset counter
  nbtot_specwin_ = 0;
  moyfc_=moytt_=0.;  
  nbpmoyttfc_=0;

  return;
}
//JEC 13/12/10 End