source: Sophya/trunk/AddOn/TAcq/branap.cc@ 4016

#include <stdlib.h>
#include <string.h>

#include "branap.h"
#include "minifits.h"
#include "strutilxx.h"
#include "sopnamsp.h"

//--------------------------------------------------------------
// Projet BAORadio - (C) LAL/IRFU  2008-2010
// Classe de gestion des parametres programmes d'analyse 
//--------------------------------------------------------------
/*!
   \class BRAnaParam
   \ingroup TAcq

   \brief parameter set for analysis programs (specmfib.cc , vismfib.cc)
*/ 

/* --Methode-- */
BRAnaParam::BRAnaParam(uint_4 nmean, uint_4 nzon, uint_4 npaqz)
{
  outpath_="./";
  fgfitsout_=false;
  imin_=imax_=0;
  istep_=1;
  rdsamefc_=true;   // read paquets with same frame counter
  freqmin_=0; freqmax_=-1;
  nbinfreq_=1;  
  paqsize_=16424;
  nzones_=nzon;
  npaqinzone_=npaqz;
  fgdatafft_=false;   fgforceraworfft_=false; 
  fgsinglechannel_=false; fgforcesingleortwochan_=false;
  prtlevel_=0;
  prtmodulo_=50000;
  prtmodulo2_=10;

  nmean_=nmean;
  fgtimeinterval_=false;
  timeinterval_=1.;
  nbloc_=1;

  nbcalgrp_=1;
  nthreads_=1;

  spec_win_sz_=0;
  spw_ext_sz_=0;
  nbmax_specwfiles_=0;

  vmin_=0.;  vmax_=9e99;
  nbands_=0; bandfirst_ = bandlast_ = 0;
  fgdtpaq_ = fgdtms_ = false;

  fgfreqfilter_=false; //JEC 1/2/11
  medhalfwidth_=50;    //JEC 6/4/11

  gainfile_="";

  proctimeduration_=9.e9;
  fgtimeselect_=false;
}

/* --Methode-- */
int BRAnaParam::DecodeArgs(int narg, char* arg[])
{
  if ((narg>1)&&(strcmp(arg[1],"-h")==0))  return Usage(false);
  if (narg<5) return Usage(true);

  bool okarg=false;
  int ka=1;
  while (ka<(narg-1)) {
    if (strcmp(arg[ka],"-act")==0) {
      action_=arg[ka+1];
      ka+=2;
    }
    else if (strcmp(arg[ka],"-out")==0) {
      outpath_=arg[ka+1];
      size_t lenp=outpath_.size();
      if ((lenp>0)&&(outpath_[lenp-1]!='/'))  outpath_+='/';
      ka+=2;
    }
    else if (strcmp(arg[ka],"-fitsout")==0) {
      fgfitsout_=true; 
      ka++;
    }
    else if (strcmp(arg[ka],"-nmean")==0) {
      nmean_=atoi(arg[ka+1]);
      ka+=2;
    }
    else if (strcmp(arg[ka],"-tmint")==0) {
      timeinterval_=atof(arg[ka+1]);
      fgtimeinterval_=true;  ka+=2;
    }
    else if (strcmp(arg[ka],"-nbloc")==0) {
      nbloc_=atoi(arg[ka+1]);
      ka+=2;
    }
    else if (strcmp(arg[ka],"-freq")==0) {
      sscanf(arg[ka+1],"%d,%d,%d",&freqmin_,&freqmax_,&nbinfreq_);  
      ka+=2;
    }
    else if (strcmp(arg[ka],"-zones")==0) {
      int nzon=4;
      int npaqz=128;
      sscanf(arg[ka+1],"%d,%d",&nzon,&npaqz);
      nzones_=nzon;  npaqinzone_=npaqz;
      ka+=2;
    }
    else if (strcmp(arg[ka],"-prt")==0) {
      sscanf(arg[ka+1],"%d,%ld,%ld",&prtlevel_,&prtmodulo_,&prtmodulo2_);
      ka+=2;
    }
    else if (strcmp(arg[ka],"-nthr")==0) {
      nthreads_=atoi(arg[ka+1]);
      ka+=2;
    }
    else if (strcmp(arg[ka],"-nvcal")==0) {
      nbcalgrp_=atoi(arg[ka+1]);
      ka+=2;
    }
    else if (strcmp(arg[ka],"-nosfc")==0) {
      rdsamefc_=false;
      ka++;
    }
    else if (strcmp(arg[ka],"-singlechan")==0) {
      fgsinglechannel_=true; fgforcesingleortwochan_ = true; 
      ka++;
    }
    else if (strcmp(arg[ka],"-twochan")==0) {
      fgsinglechannel_=false; fgforcesingleortwochan_ = true; 
      ka++;
    }
    else if (strcmp(arg[ka],"-rawdata")==0) {
      fgforceraworfft_=true; fgdatafft_ = false; 
      ka++;
    }
    else if (strcmp(arg[ka],"-fftdata")==0) {
      fgforceraworfft_=true; fgdatafft_ = true; 
      ka++;
    }
    else if (strcmp(arg[ka],"-varcut")==0) {
      sscanf(arg[ka+1],"%lg,%lg",&vmin_,&vmax_);
      ka+=2;
    }
    else if (strcmp(arg[ka],"-nband")==0) {
      sscanf(arg[ka+1],"%d,%d,%d",&nbands_,&bandfirst_,&bandlast_);
      ka+=2;
    }
    else if (strcmp(arg[ka],"-fdtpaq")==0) {
      fgdtpaq_=true;
      ka++;
    }
    else if (strcmp(arg[ka],"-fdtms")==0) {
      fgdtms_=true;
      ka++;
    }
    else if (strcmp(arg[ka],"-freqfilter")==0) {
      fgfreqfilter_=true;  ka++;
      if (strcmp(arg[ka],"-")!=0) medhalfwidth_=atof(arg[ka]); 
      ka++;
    }    
    //-tmproc hh:mm:ss,nseconds
    else if (strcmp(arg[ka],"-tmproc")==0) {
      int ah=0,am=0;
      double as=0, als=0;
      sscanf(arg[ka+1],"%d:%d:%lg,%lg",&ah,&am,&as,&als);
      fgtimeselect_=true;  proctimeduration_=als;  
      proctimestart_.SetHour(ah,am,as); 
      ka+=2;
    }
    else if (strcmp(arg[ka],"-gain")==0) {
      gainfile_=arg[ka+1];
      ka+=2;
    }
    else if (strcmp(arg[ka],"-tspwin")==0) {
      int ai1=0,ai2=0,ai3=0;
      sscanf(arg[ka+1],"%d,%d,%d",&ai1,&ai2,&ai3);
      spec_win_sz_=ai1;  spw_ext_sz_=ai2;   nbmax_specwfiles_=ai3; 
      ka+=2;
    }
    else if (strcmp(arg[ka],"-in")==0) {
      if ((narg-ka)<4)  {
        cout << " BRAnaParam::DecodeArgs() / Argument error " << endl;
        return Usage(true);
      }
      sscanf(arg[ka+1],"%d,%d,%d",&imin_,&imax_,&istep_);  ka+=2;
      while(ka<(narg-1)) {
        string inpath = arg[ka]; 
        size_t lenp=inpath.size();
        if (lenp<1)  inpath="./";
        if ((lenp>0)&&(inpath[lenp-1]!='/'))  inpath+='/';
        vector<string> fiblist;
        string sa1 = arg[ka+1]; 
        FillVStringFrString(sa1, fiblist, ',');
        char dbuff[32];
        for(size_t i=0; i<fiblist.size(); i++) {
          sprintf(dbuff,"Fiber%d/",(int)atoi(fiblist[i].c_str()));
          dirlist_.push_back(inpath+dbuff);
        }
        ka += 2;
      }
      okarg=true;
    }
    else ka++;
  }

  if (!okarg) {
    cout << " BRAnaParam::DecodeArgs() / Argument error " << endl;
    return Usage(true);
  }
  return 0;
}

/* --Methode-- */
int BRAnaParam::Usage(bool fgshort)
{
  cout << " --- BRAnaParam : Reading/Processing BAORadio FITS files parameters " << endl;
  cout << " Usage:  prgname [-act ACT] [-out OutPath] [-fitsout] \n" 
       << "                 [-nmean NMean] [-tmint dtime] [-zones NZones,nPaqinZone] \n" 
       << "                 [-nbloc NBloc] [-freq NumFreqMin,NumFreqMax,NBinFreq] \n" 
       << "                 [-prt lev,modulo,mod2] [-nvcal n] [-nthr n] [-nosfc]\n"
       << "                 [-singlechan] [-twochan] [-fftdata] [-rawdata] \n"
       << "                 [-freqfilter medhw] [-gain filename] [-varcut min,max] [-nband nband,first,last] \n"
       << "                 [-freqfilter] [-gain filename] [-varcut min,max] [-nband nband,first,last] \n"
       << "                 [-tmproc hh:mm:ss,nseconds] [-fdtpaq] [-fdtms] [-tspwin wsz,extsz,nfiles] \n" 
       << "         -in Imin,Imax,Istep InPath FiberList [InPath2 FiberList2 InPath3 FiberList3 ...] \n" << endl; 
  if (fgshort) {
    cout << " prgname -h for detailed instructions" << endl;
    return 5;
  }
  cout << " -act Action: cube3d or vis or viscktt or mspec \n"
       << "    cube3d: create 3D fits cubes \n" 
       << "    vis: compute visibilites (vismfib program) \n" 
       << "    viscktt: compute visibilities and check TimeTag/FrameCounter (vismfib program)\n"
       << "    mspec: compute and save mean spectra for each channel \n"
       << "    bproc: run BRBaseProcessor for debug/printing (use -prt)  \n"
       << " -out OutPath: Output directory name \n"
       << " -fitsout : Force FITS format for output files  \n"
       << " -nmean NMean: Number of packet used for spectra/visibility computation \n"
       << " -tmint dtime: Time interval definition (instead of nmean) for spectra/visibility computation \n"
       << " -nbloc NBloc: Number of MemMgr blocs in output file\n"
       << " -zones NZones,NbPaqinZone : Number of Zones and number of paquets in one zone (RAcqMemZoneMgr) \n" 
       << " -freq NumFreqMin,NumFreqMax,NBinFreq : Frequency zone and number of bins \n" 
       << " -prt lev,modulo,mod2 : Print level (0,1,2...), counter modulo (10000,50000...), mod2 (10,50...) for filenums \n" 
       << " -nvcal n : number of BRVisibilityCalculator objects running in parallel in BRVisCalcGroup  (default=1) \n" 
       << " -nthr n : number of threads for parallel execution in BRVisibilityCalculator (default=1) \n" 
       << " -nosfc : Don't force reading with SAME FrameCounter \n" 
       << " -singlechan : Force one channel per fiber \n"
       << " -twochan : Force two channels per fiber \n"
       << " -rawdata : Force raw data mode (firmware raw) \n"
       << " -fftdata : Force FFT data mode (firmware fft) \n"
       << " -varcut min,max : min-max cut on variance \n" 
       << " -nband nband,first,last: numbers of freq. bands and first and last bands used for cuts\n"
       << " -tmproc hh:mm:ss,nseconds : processing time window definition \n"
       << " -fdtpaq : force per paquet data table filling (specmfib) \n"
       << " -fdtms : force time averaged/filtered data table filling; use -freq for defining bands (specmfib)\n"
       << " -freqfilter medhw: force median filtering on the frequencies \n"
       << "   with half window width medhw (use - for default=" << medhalfwidth_ << ") \n"
       << " -gain filename : spectral response fits file name \n" 
       << " -tspwin wsz,extsz,nfiles : spectra time (paquet no) window (ex: -tspwin 120,4,5) \n"
       << "    wsz,extsz= window,extension size;  nfiles= maximum number of windows saved \n"
       << " -in : input files/directory definition : \n"
       << "       Imin,Imax,Istep: fits files signalII.fits Imin<=II<=Imax Istep=increment \n" 
       << "       InPath: Input directerory fits files in InPath/FiberJJ directory \n"
       << "       FiberList: List of fiber numbers (JJ - Ex: 2,1,4 )  \n" << endl;
  return 1;
}

/* --Fonction-- */
int BRAnaParam::PaqSizeFromFits()
{
  uint_4 paqsz, npaqf;
  char flnm[1024];
  sprintf(flnm,"%s/signal%d.fits",dirlist_[0].c_str(),imin_);
  bool fgdatafft_in_fits=false;
  bool fgsinglechan_in_fits=false;
  SOPHYA::TimeStamp tmstart;
  int rc = DecodeMiniFitsHeader(flnm,paqsize_, npaqinfile_,fgdatafft_in_fits, fgsinglechan_in_fits,tmstart); 
  if (!fgforcesingleortwochan_) fgsinglechannel_=fgsinglechan_in_fits;
  if (!fgforceraworfft_) fgdatafft_=fgdatafft_in_fits;
  if (fgtimeselect_) {
    int year,month,day; 
    tmstart.GetDate(year,month,day);
    proctimestart_.SetDate(year,month,day);
    proctimeend_.Set(proctimestart_.ToDays()+proctimeduration_/86400.);
  }
  return rc;
}

/* --Fonction-- */
ostream& BRAnaParam::Print(ostream& os)
{
  os << " BRAnaParam::Print() dirlist_.size()=" << dirlist_.size() << " Input directories: " << endl;
  for(size_t k=0; k< dirlist_.size(); k++) 
    cout << k+1 << " :  " << dirlist_[k] << endl;
  cout << " IMin= " << imin_ << " IMax= " << imax_ << " IStep= " << istep_ 
       << ((rdsamefc_)?" SameFrameCounter read mode":" AllOKPaquets read mode ") << endl;
  if (fgtimeselect_) {
     cout << " Processing time window, StartTime=" << proctimestart_ << " duration= " << proctimeduration_ 
              << " EndTime=" << proctimeend_ << endl;
  }
  cout << " OutPath= " << outpath_ << (fgfitsout_?" force FITS output":" PPF output") << endl;
  cout << " Action=" << action_ << "  NMean=" << nmean_ << " NBloc=" << nbloc_ << endl;
  cout << ((fgtimeinterval_)?" Time Interval mode for visi/spectra ":" NPaq Interval mode") 
       << " TimeInterval=" << timeinterval_ << " seconds " << endl;
  cout << " FreqMin= " << freqmin_ << " FreqMax= " << freqmax_ << " NBinFreq= " << nbinfreq_ << endl;
  if (fgdtpaq_) cout<<  " Fill Per Paquet DadaTable "; 
  if (fgdtms_) cout<<  " Fill Time Averaged/Filtered binned power DataTable "; 
  if (!fgdtpaq_&&!fgdtms_)  cout << " NO DataTable " << endl;
  else cout << endl;
  cout << " GainFileName=" << gainfile_ << " Variance: Min= " << vmin_ << " Max= " << vmax_ 
       << " Bands: N=" <<  nbands_ << " First=" << bandfirst_ << " Last=" << bandlast_
       << endl;
  if (fgfreqfilter_) {
    cout << " force frequence median filtering (action gain), half width =" << medhalfwidth_
         << endl;
  } else {
    cout<<" NO freq. filtering" << endl;
  }
  cout << " Spectra TimeWindow (Nb.Paquets) Size=" << spec_win_sz_ << " ExtensionSize=" << spw_ext_sz_ 
       << " MaxNbFile=" << nbmax_specwfiles_ << endl; 
  cout << " PaqSize=" << paqsize_ << "  -  NZones=" << nzones_ << " NPaqZone=" << npaqinzone_ 
       << " PrtLevel=" << prtlevel_ << " PrtCntModulo=" <<  prtmodulo_ << "," << prtmodulo2_ << endl;
  cout << " AcqDataMode: " << ((fgdatafft_)?" Data_FFT " : " Data_Raw " ) 
       << ((fgsinglechannel_)?" SingleChannel " : " TwoChannels " ) << endl;
  cout << " NbVisibCalculator in Group: " <<  nbcalgrp_ << " with N//threads: " << nthreads_ << endl;


  return os;
}

/* --Fonction-- */
int BRAnaParam::DecodeMiniFitsHeader(const char* filename, uint_4& paqsz, uint_4& npaq, 
                                     bool& fgdatafft, bool& fgsinglechannel, SOPHYA::TimeStamp& tmstart) 
{
  cout << " DecodeMiniFitsHeader - Opening file: " << filename << endl;
  MiniFITSFile mff(filename, MF_Read);
  string acqmode=mff.GetKeyValue("ACQMODE");
  cout << "DecodeMiniFitsHeader()... Type=" << mff.DataTypeToString() << " NAxis1=" << mff.NAxis1() 
       << " NAxis2=" << mff.NAxis2() << " AcqMode=" << acqmode << endl;
  paqsz = mff.NAxis1();
  npaq = mff.NAxis2();
  if (acqmode.substr(0,3)=="fft") fgdatafft=true;
  if (acqmode.substr(3,3)=="fft") fgdatafft=true;  // pour ordfft1c ordfft2c - Cedric Viou, 15 Juin 2011
  if (acqmode.find("1c") < acqmode.length()) fgsinglechannel=true;
  string fkvs=mff.GetKeyValue("DATEOBS");
  if (fkvs.length()>0)  tmstart.Set(fkvs);
  fkvs=mff.GetKeyValue("TMSTART");
  if (fkvs.length()>0) tmstart.Set(fkvs);
  return 0;
}