//---------------------------------------------------------------- // Projet BAORadio - (C) LAL/IRFU 2008-2010 // Classes de threads de traitement pour BAORadio //---------------------------------------------------------------- #include #include #include #include #include #include #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 //--------------------------------------------------------------------- /* --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 medianFilterPaq(spec_window_.SizeX()); //we keep freq. = 0 for compatibility for (sa_size_t f=0; f val1; TVector 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 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= fMin) ? f-hww : fMin; sa_size_t last = (f+hww <= fMax) ? f+hww : fMax; std::vector 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_ 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