Changeset 3993 in Sophya

Timestamp:

May 13, 2011, 12:37:14 AM (14 years ago)

Author:

ansari

Message:

Ajout possibilite de faire un DataTable (NTuple) avec la puissance spectrale dans plusieurs bandes en fonction du temps, Reza 13/05/2011

Location:

trunk/AddOn/TAcq

Files:

• branap.cc (modified) (6 diffs)
• branap.h (modified) (2 diffs)
• brbaseproc.cc (modified) (2 diffs)
• brbaseproc.h (modified) (2 diffs)
• brproc.cc (modified) (11 diffs)
• brproc.h (modified) (4 diffs)
• brprocGain.cc (modified) (2 diffs)
• brviscalc.cc (modified) (1 diff)
• specmfib.cc (modified) (2 diffs)

trunk/AddOn/TAcq/branap.cc

@@ -22 +22 @@
   istep_=1;
   rdsamefc_=true;   // read paquets with same frame counter
-  freqmin_=freqmax_=0;
+  freqmin_=0; freqmax_=-1;
   nbinfreq_=1;  
   paqsize_=16424;
@@ -40 +40 @@
   vmin_=0.;  vmax_=9e99;
   nbands_=0; bandfirst_ = bandlast_ = 0;
-  fgdt_ = false;
+  fgdtpaq_ = fgdtms_ = false;
 
   fgfreqfilter_=false; //JEC 1/2/11
@@ -133 +133 @@
       ka+=2;
     }
-    else if (strcmp(arg[ka],"-filldt")==0) {
-      fgdt_=true;
+    else if (strcmp(arg[ka],"-fdtpaq")==0) {
+      fgdtpaq_=true;
+      ka++;
+    }
+    else if (strcmp(arg[ka],"-fdtms")==0) {
+      fgdtms_=true;
       ka++;
     }
@@ -205 +209 @@
        << "                 [-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] [-filldt] [-tspwin wsz,extsz,nfiles] \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) {
@@ -232 +236 @@
        << " -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"
+       << " -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"
-       << " -filldt : force data table filling \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"
@@ -283 +288 @@
   cout << " Action=" << action_ << "  NMean=" << nmean_ << " NBloc=" << nbloc_ << endl;
   cout << " FreqMin= " << freqmin_ << " FreqMax= " << freqmax_ << " NBinFreq= " << nbinfreq_ << endl;
-  cout<<  ((fgdt_)?" Fill DadaTable ":" NO DataTable") << 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_

trunk/AddOn/TAcq/branap.h

@@ -58 +58 @@
   double vmin_,vmax_;   // coupure min-max sur la variance 
   int nbands_,bandfirst_,bandlast_; //bandes en frequences pour mean, variance et cut
-  bool fgdt_; //if true fill datatable
+  bool fgdtpaq_; //if true fill per per paquet datatable
+  bool fgdtms_; //if true fill time avergared/filtered spectral power  datatable
 
   bool fgfreqfilter_; // if true force median filtering on the frequencies
@@ -65 +66 @@
   string gainfile_;     // nom du fichier de gain (reponse spectrale)
 
+  // Remplissage DataTable power = f(time) 
+  
   // definition fenetre en temps pour le traitement des paquets 
   SOPHYA::TimeStamp proctimestart_; //  temps de debut de la fenetre

trunk/AddOn/TAcq/brbaseproc.cc

@@ -40 +40 @@
   fgokallfibers_=true;
   totprocnpaq_=0;
+  startdate_=0;
   setNameId("baseproc",0); 
   ClearProcTimeWindow();
@@ -119 +120 @@
       }
       cts_=memgr_.GetAuxData(mid_)->FillTime();   // get associated date/time (DATEOBS)
+      if (startdate_==0)  startdate_=cts_.DaysPart();
+
       for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
         fbuff_[fib] = memgr_.GetMemZone(mid_,fib);

trunk/AddOn/TAcq/brbaseproc.h

@@ -50 +50 @@
   // get Observation Time ( obtained from MemZoneMgr )
   inline SOPHYA::TimeStamp& getObsTime()  { return cts_; }
+  // get Obsservation time in seconds since the start date 
+  // Temps ecoule en seconde a partir de 0h00 de la date du premier fichier
+  inline double getObsTimeSeconds() { return ((cts_.DaysPart()-startdate_)*86400.+cts_.SecondsPart()); }
+
   // Attention : les methodes suivantes ne sont pas protegees pour le numero de fibre
   // Renvoie le numero de frame-counter courant - si fgz=true, soustrait le FC du premier paquet 
@@ -89 +93 @@
   uint_8 totprocnpaq_; 
   SOPHYA::TimeStamp cts_;   // current time stamp, get from MemZoneMgr
+  int_8 startdate_;         // Temps de depart cts_.DaysPart() 
 
   string bpnom_;    // nom du processeur

trunk/AddOn/TAcq/brproc.cc

@@ -73 +73 @@
   numfile_=0;
   totnbpaq_=0;
+  moyfc_=moytt_=0.;
+  nbpmoyttfc_=0;
 
   size_t nchan=(fgsinglechannel_?memgr_.NbFibres():2*memgr_.NbFibres());
@@ -79 +81 @@
     nbpaq4mean_[i]=nbadpaq_[i]=0;
   }
-
-  //
   
 
@@ -92 +92 @@
   ofsdtp_=NULL;
   dtp_=NULL;
-  
-  //  cout << "(JEC) creation tuple de " << nchan*(1+numberOfBands_) << " doubles " << endl;
-  //  xnt_=new double[nchan*(1+numberOfBands_)];  // CHECK : ATTENTION la taille depend de nombre de colonne du NTuple !
   xnt_=NULL;
+  ofsdtpms_=NULL;
+  dtpms_=NULL;
+  xntms_=NULL;
+  firstfreqms_=0; 
+  lastfreqms_=-1;
+  nbandms_=-1;
 }
 
@@ -113 +116 @@
     delete ofsdtp_;
   }
+  if (dtpms_) {
+    cout << *dtpms_; 
+    delete dtpms_;
+    delete ofsdtpms_;
+  }
   if (xnt_)   delete xnt_;
+  if (xntms_)   delete xntms_;
   cout << " ------------------------------------------------------------------------ " << endl;
 }
@@ -142 +151 @@
 
 /* --Methode-- */
-void BRMeanSpecCalculator::DefineDataTable()
-{
-  cout << "(JEC) BRMeanSpecCalculator::DefineDataTable START" << endl;
-  string dtfile="!"+outpath_+"/dtspec.fits";
+void BRMeanSpecCalculator::DefinePerPaquetDataTable()
+{
+  cout << "BRMeanSpecCalculator::DefinePerPaquetDataTable: Creating DataTable" << endl;
+  string dtfile="!"+outpath_+"/dtspecpaq.fits";
   ofsdtp_ = new FitsInOutFile(dtfile,FitsInOutFile::Fits_Create);
   dtp_ = new SwFitsDataTable(*ofsdtp_,1024,true);
   char cnom[32];
-  size_t nchan=(fgsinglechannel_?memgr_.NbFibres():2*memgr_.NbFibres());
+  size_t nchan=mspecmtx_.NRows();
   for(int i=0; i<nchan; i++) {
     sprintf(cnom,"variance%d",i);
@@ -161 +170 @@
   }
   /*
-  for(int i=0; i<nchan; i++) {
+    for(int i=0; i<nchan; i++) {
     sprintf(cnom,"sigma%d",i);
     dtp_->AddFloatColumn(cnom);    
-  }
+    }
   */
   //  xnt_=new double[nchan*2];  CHECK : faut-il reallouer ?
-  cout << "(JEC) creation tuple de " << nchan*(1+numberOfBands_) << " doubles " << endl;
+  cout << " ...Number of columns " << nchan*(1+numberOfBands_) << " doubles -> to FITS file" << dtfile << endl;
   xnt_=new double[nchan*(1+numberOfBands_)];
-  cout << "(JEC) BRMeanSpecCalculator::DefineDataTable END" << endl;
+
+}
+
+/* --Methode-- */
+void BRMeanSpecCalculator::DefineTimeAvgPowerDataTable(int first, int last, int nband)
+{
+  cout << "BRMeanSpecCalculator::DefineTimeAvgPowerDataTable: Creating time averged power spectrum" << endl;
+  string dtfile="!"+outpath_+"/avgspec.fits";
+  ofsdtpms_ = new FitsInOutFile(dtfile,FitsInOutFile::Fits_Create);
+  dtpms_ = new SwFitsDataTable(*ofsdtpms_,1024,true);
+  dtpms_->AddDoubleColumn("mtimetag");
+  dtpms_->AddDoubleColumn("time");
+  
+  if (last>first) {
+    if ((first<0)||(first>=mspecmtx_.NCols()))  first=0;
+    if ((last<0)||(last>mspecmtx_.NCols()))  last=mspecmtx_.NCols();
+    if (nband<1)  nband=1;
+  }
+  else {
+    first=1;  last=mspecmtx_.NCols();  nband=1;
+  }
+  firstfreqms_=first;   lastfreqms_=last;  nbandms_=nband;
+  char cnom[32];
+  size_t nchan=mspecmtx_.NRows();
+  for(int i=0; i<nchan; i++) {
+    for(int j=0; j<nband; j++) {
+      sprintf(cnom,"pow%dc%d",j,i);
+      dtpms_->AddFloatColumn(cnom); 
+    }   
+  }
+  cout << " ...Number of columns " << (2+nchan*nband) << " doubles -> to FITS file" << dtfile << endl;
+  xntms_=new double[nchan*nband+2];
 }
   
@@ -181 +221 @@
   ibandfirst_=ibandfirst; ibandlast_=ibandlast;
 
-  cout << "(JEC): SetNumberOfBands (END) : "
+  cout << "BRMeanSpecCalculator::SetNumberOfBands (END) : "
        << numberOfBands_ << " "
        << ibandfirst_ << " "
@@ -264 +304 @@
   }
 
+  // Pour calcul MeanTimeTag , MeanFrameCounter 
+  moyfc_+=curfc_[0];
+  moytt_+=((double)vpaq_[0].TimeTag()/1.25e8);
+  nbpmoyttfc_++;
   totnbpaq_++;
   return 0;
@@ -371 +415 @@
     }
   }
+  mspecmtx_.Info()["ENDTIME"]=getObsTime();
+  sigspecmtx_.Info()["ENDTIME"]=getObsTime();
+  if (nbpmoyttfc_>1)  {
+    moyfc_/=(double)nbpmoyttfc_;
+    moytt_/=(double)nbpmoyttfc_;
+  }
+  string ikey,ikdesc;    
+  ikey="MeanFC";  ikdesc="Mean FrameCounter";
+  mspecmtx_.Info().SetI(ikey,moyfc_); 
+  mspecmtx_.Info().SetComment(ikey,ikdesc);
+  sigspecmtx_.Info().SetI(ikey,moyfc_); 
+  sigspecmtx_.Info().SetComment(ikey,ikdesc);
+  
+  ikey="MeanTT";  ikdesc="Mean TimeTag";
+  mspecmtx_.Info().SetD(ikey,moytt_);
+  mspecmtx_.Info().SetComment(ikey,ikdesc);
+  sigspecmtx_.Info().SetD(ikey,moytt_);
+  sigspecmtx_.Info().SetComment(ikey,ikdesc);
+  
   char nfile[64];
   string flnm;
@@ -390 +453 @@
   numfile_++;  
 
+  if (dtpms_!=NULL)  FillPwrTmDTable(mspecmtx_);
+
   for(size_t i=0; i<nbpaq4mean_.size(); i++) nbpaq4mean_[i]=0;
   mspecmtx_ = (r_4)(0.);
   sigspecmtx_ = (r_4)(0.);
+  moyfc_=moytt_=0.;  
+  nbpmoyttfc_=0;
   return;
 }
@@ -406 +473 @@
   fos << spec_window_;
   cout << " SaveSpectraWindow() " << nbtot_specwin_ << "- file " << nfile << " created " << endl;
+}
+
+/* --Methode-- */
+void BRMeanSpecCalculator::FillPwrTmDTable(TMatrix< r_4 >& specmtx)
+{
+  if (dtpms_==NULL) return;
+  xntms_[0]=moytt_;  
+  xntms_[1]=getObsTimeSeconds();
+  size_t knt=2;  
+  sa_size_t wband=(lastfreqms_-firstfreqms_)/nbandms_;
+  for(sa_size_t ir=0; ir<specmtx.NRows(); ir++) {
+    for(sa_size_t jf=0; jf<specmtx.NCols()-wband; jf+=wband) {
+      TVector<r_4> srow=specmtx.Row(ir);
+      xntms_[knt]=srow.SubVector(Range(jf,jf+wband-1)).Sum();
+      knt++;
+    }
+  }
+  dtpms_->AddRow(xntms_);
 }
 

trunk/AddOn/TAcq/brproc.h

@@ -63 +63 @@
   inline void SetMaxNbSpecWinFiles(uint_4 nmax=0)  {  nmaxfiles_specw_=nmax; }
 
-// Pour definir le datatable a remplir - pas de DataTable rempli si pas appele 
-  virtual void DefineDataTable();
+// Pour definir le datatable (moyenne spectre paquet par paquet)  
+// Le DataTable est rempli par FlagBadPackets - Pas de DataTable rempli si pas appele 
+  virtual void DefinePerPaquetDataTable(); 
+
+// Pour definir le datatable puissance spectrale moyennee ou filtree en fonction du temps
+  virtual void DefineTimeAvgPowerDataTable(int first=1, int last=-1, int nband=1); 
 
   inline void SetVarianceLimits(double vmin=0., double vmax=9.e99)
@@ -85 +89 @@
   virtual void SaveMeanSpectra();  // Pour sauver les spectres moyennes ( + sigma ) dans un fichier 
   virtual void SaveSpectraWindow();  // Pour sauver les spectres de la fenetre temporel ds un fichier 
-
+  virtual void FillPwrTmDTable(TMatrix< r_4 >& specmtx);  // pour remplir la puissance spectrale en fonction du temps
    
   uint_4 nmean_;  // Nombre de spectres pour le calcul des moyennes 
@@ -95 +99 @@
   uint_8 totnbpaq_;
 
-  TMatrix< r_4 > mspecmtx_;    // Matrice des spectres moyennees 
+  TMatrix< r_4 > mspecmtx_;      // Matrice des spectres moyennees 
   TMatrix< r_4 > sigspecmtx_;    // Matrice des sigmas des spectres  
+  double moyfc_, moytt_;         // FrameCounter/TimeTag moyen pour chaque matrice 
+  uint_8 nbpmoyttfc_;            // Nombre de TimeTag/FrameCounter somme  
 
   TArray< r_4 > spec_window_;    // fenetre en temps sur les spectres  
@@ -117 +123 @@
   vector<uint_8> nbadpaq_;
 
+  // DataTable for paquet/paquet mean/variance for frequency bands 
   FitsInOutFile* ofsdtp_;   // Output fits stream for datatable 
-  SwFitsDataTable* dtp_;    // DataTable 
+  SwFitsDataTable* dtp_;    // per paquet DataTable 
   double* xnt_; 
+  // DataTable for time averaged spectrum mean per frequency band
+  FitsInOutFile* ofsdtpms_;   // Output fits stream for datatable for time averaged spectrum  
+  SwFitsDataTable* dtpms_;    // DataTable for time averaged spectrum  
+  double* xntms_; 
+  int firstfreqms_, lastfreqms_, nbandms_;
 }; 
 

trunk/AddOn/TAcq/brprocGain.cc

@@ -185 +185 @@
     medfiltspecmtx_.Info()["NPAQSUM"] = nbtot_specwin_;
     medfiltspecmtx_.Info()[buff] = nbtot_specwin_;
-    //JEC timeStamp in the day @ millisec: faudra gerer le passage a minuit!
-    medfiltspecmtx_.Info()["TIMEWIN"] = (windowTime_.SecondsPart())*1000.;
-
-
     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;
@@ -206 +217 @@
   //increment the file index
   nummedianfile_++;
+
+  if (dtpms_!=NULL)  FillPwrTmDTable(mspecmtx_);
     
   //reset the matirx
   medfiltspecmtx_ = (r_4)(0.);
-
   //reset counter
   nbtot_specwin_ = 0;
+  moyfc_=moytt_=0.;  
+  nbpmoyttfc_=0;
 
   return;

trunk/AddOn/TAcq/brviscalc.cc

@@ -165 +165 @@
 int BRVisibilityCalculator::SelectFreqBinning(uint_4 freq1, uint_4 freq2, uint_4 nbfreq)
 {
+  if (freq2<freq1) {
+    freq1=0; freq2=vismtx_.NCols()-1; nbfreq=1; 
+  }
   jf1_=freq1; jf2_=freq2;
   if ((jf1_<1)||(jf1_>=vismtx_.NCols()))  jf1_=1;

trunk/AddOn/TAcq/specmfib.cc

@@ -85 +85 @@
     //JEC 27/1/11 see below    if (par.gainfile_.length()>0) procms.ReadGainFitsFile(par.gainfile_);
     procms.SetPrintLevel(par.prtlevel_,par.prtmodulo_);
-    if (par.fgdt_) procms.DefineDataTable();
+    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_);
 
@@ -99 +101 @@
         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_);
     }