source: Sophya/trunk/AddOn/TAcq/specmfib.cc@ 4070

// Utilisation de SOPHYA pour faciliter les tests ...
#include "sopnamsp.h"
#include "machdefs.h"

/* ---------------------------------------------------------- 
   Programme de lecture multifibres, calcul de spectres 
   BAORadio -   LAL/IRFU      R. Ansari, C. Magneville, J.E. Campagne
   2010 - 2011
   ---------------------------------------------------------- */
/*!
   \ingroup TAcq
   \file specmfib.cc

   \brief Data cleaning and computation of mean spectrum of gain/normalization spectrum

   This program uses the BRMeanSpecCalculator, BRGainCalculator and BRMultiFitsReader classes
   to clean data and compute mean spectra and normalization/gain spectrum from 
   raw or FFT data read from FITS files (data in BRPaquet format). 
   The BRFFTCalculator class is used to compute fourier coefficients on raw data and 
   BRAnaParam class is used for defining and decoding the program run parameters.

  \verbatim
csh> specmfib 
 ===> specmfib.cc: decoding command line arguments 
 --- BRAnaParam : Reading/Processing BAORadio FITS files parameters 
 Usage:  prgname [-act ACT] [-out OutPath] [-fitsout] 
                 [-nmean NMean] [-tmint dtime] [-zones NZones,nPaqinZone] 
                 [-nbloc NBloc] [-freq NumFreqMin,NumFreqMax,NBinFreq] 
                 [-prt lev,modulo,mod2] [-nvcal n] [-nthr n] [-nosfc]
                 [-singlechan] [-twochan] [-fftdata] [-rawdata] 
                 [-freqfilter medhw] [-gain filename] [-varcut min,max] [-nband nband,first,last] 
                 [-freqfilter] [-gain filename] [-varcut min,max] [-nband nband,first,last] 
                 [-tmproc hh:mm:ss,nseconds] [-fdtpaq] [-fdtms] [-tspwin wsz,extsz,nfiles] 
         -in Imin,Imax,Istep InPath FiberList [InPath2 FiberList2 InPath3 FiberList3 ...] 

 prgname -h for detailed instructions
  \endverbatim 
*/

// include standard c/c++
#include <iostream>
#include <string>
#include <exception>

// include sophya mesure ressource CPU/memoire ...
#include "resusage.h"
#include "ctimer.h"
#include "timing.h"
#include "timestamp.h"
#include "strutilxx.h"
#include "tarrinit.h"
#include "histinit.h"

#include "brpaqu.h"
#include "brfitsrd.h"
#include "brproc.h"

//JEC 19/1/11 START
#include "brprocGain.h"
//JEC END

#include "brdiskw.h"

#include "branap.h"



//----------------------------------------------------
//----------------------------------------------------
int main(int narg, char* arg[])
{

  TArrayInitiator  _inia;

  int rc = 0;
  try {
    // Decodage parametres 
    BRAnaParam par;
    par.action_="mspec";
    cout << " ===> specmfib.cc: decoding command line arguments " << endl;
    rc = par.DecodeArgs(narg, arg);
    if (rc) return rc;
    rc = par.PaqSizeFromFits();
    if (rc) return rc;
    par.Print(cout);
    if ((par.action_!="cube3d")
        &&(par.action_!="mspec")
        &&(par.action_!="gain") //JEC 19/1/11 add gain action
        &&(par.action_!="bproc") 
        ) {
      cout << " !!! specmfib.cc  BAD action = " << par.action_ << " possible values: mspec,cube3d,bproc" << endl;
      return 5;
    }   

    cout << " ---------- specmfib.cc Start - Action= " << par.action_ << " ------------- " << endl;
    ResourceUsage resu;
    BRPaquet paq(par.paqsize_);
    uint_4 procsz=sizeof(float)*(paq.DataSize()+4);
//    if ((par.fgdatafft_)||(par.action_=="cube3d")) procsz = 0;
    cout << " specmfib: Creating MemZoneMgr/threads - PaqSz= " << par.paqsize_ 
         << " ProcPaqSz=" << procsz << endl;
    RAcqMemZoneMgr mmgr(par.nzones_, par.dirlist_.size(), par.npaqinzone_, par.paqsize_, procsz);
    if (par.action_ == "cube3d")  mmgr.SetFinalizedMask((uint_4)MemZS_Saved);
    else {
      if (par.fgdatafft_)  mmgr.SetFinalizedMask((uint_4)MemZS_ProcA);
      else mmgr.SetFinalizedMask((uint_4)MemZS_ProcB);
    }
    BRMultiFitsReader reader(mmgr, par.dirlist_, par.rdsamefc_, par.imin_, par.imax_, par.istep_);
    reader.SetPrintLevel(par.prtlevel_,par.prtmodulo_);

    BRMeanSpecCalculator procms(mmgr, par.outpath_, par.nmean_, par.fgdatafft_, par.fgsinglechannel_);
    procms.SetVarianceLimits(par.vmin_, par.vmax_);
    if (par.nbands_>0) procms.SetNumberOfBands(par.nbands_,par.bandfirst_,par.bandlast_);
    //JEC 27/1/11 see below    if (par.gainfile_.length()>0) procms.ReadGainFitsFile(par.gainfile_);
    procms.SetPrintLevel(par.prtlevel_,par.prtmodulo_);
    if (par.fgdtpaq_) procms.DefinePerPaquetDataTable();
    if (par.fgdtms_) procms.DefineTimeAvgPowerDataTable(par.freqmin_,par.freqmax_,par.nbinfreq_);

    if(par.fgtimeselect_) procms.SetProcTimeWindow(par.proctimestart_,par.proctimeend_);

    BRGainCalculator* procgain_p=NULL;
    if (par.action_=="gain") {
      procgain_p = 
        new BRGainCalculator(mmgr, par.outpath_, par.nmean_, par.fgdatafft_, par.fgsinglechannel_);
      procgain_p->SetPrintLevel(par.prtlevel_,par.prtmodulo_);
      if(par.fgfreqfilter_){ 
        procgain_p->SetOnMedianFreqFilt();
        //add possibility to set the half width of the Freq. Median Filter
        //default is 50 if par.fgfreqfilter_=true
        if (par.medhalfwidth_>0) procgain_p->SetMedianFilterHalfWidth(par.medhalfwidth_);
      }
      if (par.fgdtpaq_) procgain_p->DefinePerPaquetDataTable();
      if (par.fgdtms_) procgain_p->DefineTimeAvgPowerDataTable(par.freqmin_,par.freqmax_,par.nbinfreq_);
      if(par.fgtimeselect_) procgain_p->SetProcTimeWindow(par.proctimestart_,par.proctimeend_);
    }


    BRFFTCalculator procfft(mmgr, par.fgsinglechannel_);
    if(par.fgtimeselect_) procfft.SetProcTimeWindow(par.proctimestart_,par.proctimeend_);

    if (!par.fgdatafft_) { 
      procms.SetMemZAction(MemZA_ProcB);
      if(procgain_p) procgain_p->SetMemZAction(MemZA_ProcB);
    }

    if (par.spec_win_sz_>0) { 
      procms.SetSpectraWindowSize(par.spec_win_sz_, par.spw_ext_sz_);
      procms.SetMaxNbSpecWinFiles(par.nbmax_specwfiles_);
      if(procgain_p){
        procgain_p->SetSpectraWindowSize(par.spec_win_sz_, par.spw_ext_sz_);
        procgain_p->SetMaxNbSpecWinFiles(par.nbmax_specwfiles_);
      }
    }
    //JEC 27/1/11 should be done after SetSpectraWindowSize to perform x-checking
    if (par.gainfile_.length()>0) {
      procms.ReadGainFitsFile(par.gainfile_);
      if(procgain_p){ //JEC 5/4/2011 we can load previous gain...
        procgain_p->ReadGainFitsFile(par.gainfile_);
      }
    }

    FitsCubeWriter wrt(mmgr, par.outpath_, par.nbloc_);

    BRBaseProcessor* bproc_p=NULL; 
    if (par.action_=="bproc") { 
       bproc_p = new BRBaseProcessor(mmgr);
       bproc_p->SetPrintLevel(par.prtlevel_,par.prtmodulo_);
       if (!par.fgdatafft_)  bproc_p->SetMemZAction(MemZA_ProcB);
       if(par.fgtimeselect_) bproc_p->SetProcTimeWindow(par.proctimestart_,par.proctimeend_);
    }

    cout << " specmfib: Starting threads (reader + "<< par.action_ << ")" << endl;
    reader.start();
    if (par.action_ == "cube3d")   wrt.start();
    else {  // Calcul spectre moyenne 
      if (!par.fgdatafft_)  procfft.start(); 
      //JEC 19/1/11 START add the gain process alternative
      if (par.action_=="bproc")  bproc_p->start();
      else if (par.action_=="gain") procgain_p->start();
      else procms.start();
      //JEC END
    }
    usleep(200000);
    reader.join();
    if (par.action_ == "cube3d")  wrt.join();
    else {
      if (!par.fgdatafft_)  procfft.join(); 

      //JEC 19/1/11 START add the gain process alternative
      if (par.action_=="bproc")  bproc_p->join();
      else if (par.action_=="gain") procgain_p->join();
      else procms.join();
      //JEC END
    }
    mmgr.Print(cout);
    if (bproc_p) delete bproc_p;
    if (procgain_p) delete procgain_p;
    cout << resu ; 
  }
  catch (std::exception& sex) {
    cerr << "\n specmfib.cc std::exception :"  << (string)typeid(sex).name() 
         << "\n msg= " << sex.what() << endl;
    rc = 78;
  }
  catch (...) {
    cerr << " specmfib.cc catched unknown (...) exception  " << endl; 
    rc = 79; 
  } 

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

}