source: Sophya/trunk/AddOn/TAcq/brproc.cc@ 4062

Last change on this file since 4062 was 4016, checked in by ansari, 14 years ago

Ajout de commentaires d'autodocumentation Doxygen, Reza 12/08/2011
File size: 44.0 KB
Line 
1//----------------------------------------------------------------
2// Projet BAORadio - (C) LAL/IRFU 2008-2011
3// Classes de threads de traitement pour BAORadio
4//----------------------------------------------------------------
5
6#include <stdlib.h>
7#include <string.h>
8#include <unistd.h>
9#include <fstream>
10#include <signal.h>
11
12#include "pexceptions.h"
13#include "tvector.h"
14#include "ntuple.h"
15#include "datatable.h"
16#include "histos.h"
17#include "fioarr.h"
18#include "matharr.h"
19#include "timestamp.h"
20#include "ctimer.h"
21#include "fftpserver.h"
22#include "fitsarrhand.h"
23
24#include "FFTW/fftw3.h"
25
26
27#include "pciewrap.h"
28#include "brpaqu.h"
29#include "brproc.h"
30
31
32
33//---------------------------------------------------------------------
34// Classe BRMeanSpecCalculator de traitement de spectres -
35// Calcul de spectres moyennes,variance / voie + nettoyage
36// Implementation de traitement par fenetres temps-frequence
37// Le temps correspond au numero de paquet
38//---------------------------------------------------------------------
39/*!
40 \class BRMeanSpecCalculator
41 \ingroup TAcq
42
43 \brief Computes mean spectra on FFT data.
44*/
45
46/* --Methode-- */
47BRMeanSpecCalculator::BRMeanSpecCalculator(RAcqMemZoneMgr& memgr, string outpath, uint_4 nmean,
48 bool fgdatafft, bool fgsinglechan)
49 : BRBaseProcessor(memgr), outpath_(outpath), nmean_(nmean),
50 fgdatafft_(fgdatafft), fgsinglechannel_(fgsinglechan),
51 nbpaq4mean_(fgsinglechan?memgr_.NbFibres():2*memgr_.NbFibres()),
52 nbadpaq_(fgsinglechan?memgr_.NbFibres():2*memgr_.NbFibres())
53{
54 setNameId("meanSpecCalc",1);
55
56 uint_4 nb_octets_entrop = 0; //this value is valid for Dec. 2010 data at Nancay
57 const char* venvp = NULL;
58 venvp=getenv("BRANA_NBYTECUT");
59 if (venvp!=NULL){
60 nb_octets_entrop = atoi(venvp);
61 cout << "BRMeanSpecCalculator : BRANA_NBYTECUT : " << nb_octets_entrop << endl;
62 }
63
64 BRPaquet paq(memgr_.PaqSize()-nb_octets_entrop);
65
66 if (fgsinglechannel_) {
67 mspecmtx_.SetSize(memgr_.NbFibres(), paq.DataSize()/2);
68 sigspecmtx_.SetSize(memgr_.NbFibres(), paq.DataSize()/2);
69 sgain_.SetSize(memgr_.NbFibres(), paq.DataSize()/2);
70 }
71 else {
72 mspecmtx_.SetSize(2*memgr_.NbFibres(), paq.DataSize()/4);
73 sigspecmtx_.SetSize(2*memgr_.NbFibres(), paq.DataSize()/4);
74 sgain_.SetSize(2*memgr_.NbFibres(), paq.DataSize()/4);
75 }
76 mspecmtx_=(r_4)(0.);
77 sigspecmtx_=(r_4)(0.);
78 sgain_=(r_4)(1.); // Gain en fonction de la frequence, à 1 par defaut
79
80 numfile_=0;
81 totnbpaq_=0;
82 moyfc_=moytt_=0.;
83 nbpmoyttfc_=0;
84
85 size_t nchan=(fgsinglechannel_?memgr_.NbFibres():2*memgr_.NbFibres());
86
87 for(size_t i=0; i<nchan; i++) {
88 nbpaq4mean_[i]=nbadpaq_[i]=0;
89 }
90
91
92 // Definition des tailles de fenetres de spectres, etc ...
93 SetSpectraWindowSize();
94 SetMaxNbSpecWinFiles();
95 nbtot_specwin_=0;
96 SetVarianceLimits();
97 SetNumberOfBands();
98
99 ofsdtp_=NULL;
100 dtp_=NULL;
101 xnt_=NULL;
102 ofsdtpms_=NULL;
103 dtpms_=NULL;
104 xntms_=NULL;
105 firstfreqms_=0;
106 lastfreqms_=-1;
107 nbandms_=-1;
108}
109
110/* --Methode-- */
111BRMeanSpecCalculator::~BRMeanSpecCalculator()
112{
113 cout << " ---------------- BRMeanSpecCalculator()_Finalizing -------------------- " << endl;
114 uint_8 npqm=0;
115 for(size_t i=0; i<nbpaq4mean_.size(); i++) npqm+=nbpaq4mean_[i];
116 if (npqm>nmean_*nbpaq4mean_.size()/10) SaveMeanSpectra();
117 for(size_t i=0; i<nbadpaq_.size(); i++) {
118 cout << " Channel " << i << " NBadPaq=" << nbadpaq_[i] << " / TotNbPaq=" << totnbpaq_ << endl;
119 }
120 if (dtp_) {
121 cout << *dtp_;
122 delete dtp_;
123 delete ofsdtp_;
124 }
125 if (dtpms_) {
126 cout << *dtpms_;
127 delete dtpms_;
128 delete ofsdtpms_;
129 }
130 if (xnt_) delete xnt_;
131 if (xntms_) delete xntms_;
132 cout << " ------------------------------------------------------------------------ " << endl;
133}
134
135/* --Methode-- */
136void BRMeanSpecCalculator::SetSpectraWindowSize(uint_4 winsz, uint_4 wszext)
137{
138 if (winsz < 3) {
139 winsz=1; wszext=0;
140 }
141 if (wszext>=winsz/2) wszext=winsz/2;
142 sa_size_t sz[5]={0,0,0,0,0};
143 sz[0]=mspecmtx_.NCols();
144 sz[1]=mspecmtx_.NRows();
145 sz[2]=winsz+2*wszext;
146 spec_window_.SetSize(3,sz);
147 spwin_ext_sz_=wszext;
148 sz[0]=mspecmtx_.NRows();
149 sz[1]=winsz+2*wszext;
150 clnflg_.SetSize(2,sz);
151 cout << "BRMeanSpecCalculator::SetSpectraWindowSize()/Info: SpectraWindowSize()=" << GetSpectraWindowSize()
152 << " ExtensionSize=" << GetSpecWinExtensionSize() << " Overlap=" << GetSpecWinOverlapSize()
153 << " ArraySize=" << spec_window_.SizeZ() << endl;
154
155 paqnum_w_start=spwin_ext_sz_; // premiere initialisation du numero de paquet
156 return;
157}
158
159/* --Methode-- */
160void BRMeanSpecCalculator::DefinePerPaquetDataTable()
161{
162 cout << "BRMeanSpecCalculator::DefinePerPaquetDataTable: Creating DataTable" << endl;
163 string dtfile="!"+outpath_+"/dtspecpaq.fits";
164 ofsdtp_ = new FitsInOutFile(dtfile,FitsInOutFile::Fits_Create);
165 dtp_ = new SwFitsDataTable(*ofsdtp_,1024,true);
166 char cnom[32];
167 size_t nchan=mspecmtx_.NRows();
168 for(int i=0; i<nchan; i++) {
169 sprintf(cnom,"variance%d",i);
170 dtp_->AddFloatColumn(cnom);
171 }
172 for(int i=0; i<nchan; i++) {
173 for(int j=0;j<numberOfBands_;j++){
174 sprintf(cnom,"varnormb%d%d",i,j);
175 dtp_->AddFloatColumn(cnom);
176 }
177 }
178 /*
179 for(int i=0; i<nchan; i++) {
180 sprintf(cnom,"sigma%d",i);
181 dtp_->AddFloatColumn(cnom);
182 }
183 */
184 // xnt_=new double[nchan*2]; CHECK : faut-il reallouer ?
185 cout << " ...Number of columns " << nchan*(1+numberOfBands_) << " doubles -> to FITS file" << dtfile << endl;
186 xnt_=new double[nchan*(1+numberOfBands_)];
187
188}
189
190/* --Methode-- */
191void BRMeanSpecCalculator::DefineTimeAvgPowerDataTable(int first, int last, int nband)
192{
193 cout << "BRMeanSpecCalculator::DefineTimeAvgPowerDataTable: Creating time averged power spectrum" << endl;
194 string dtfile="!"+outpath_+"/avgspec.fits";
195 ofsdtpms_ = new FitsInOutFile(dtfile,FitsInOutFile::Fits_Create);
196 dtpms_ = new SwFitsDataTable(*ofsdtpms_,1024,true);
197 dtpms_->AddDoubleColumn("mtimetag");
198 dtpms_->AddDoubleColumn("time");
199
200 if (last>first) {
201 if ((first<0)||(first>=mspecmtx_.NCols())) first=0;
202 if ((last<0)||(last>mspecmtx_.NCols())) last=mspecmtx_.NCols();
203 if (nband<1) nband=1;
204 }
205 else {
206 first=1; last=mspecmtx_.NCols(); nband=1;
207 }
208 firstfreqms_=first; lastfreqms_=last; nbandms_=nband;
209 char cnom[32];
210 size_t nchan=mspecmtx_.NRows();
211 for(int i=0; i<nchan; i++) {
212 for(int j=0; j<nband; j++) {
213 sprintf(cnom,"pow%dc%d",j,i);
214 dtpms_->AddFloatColumn(cnom);
215 }
216 }
217 cout << " ...Number of columns " << (2+nchan*nband) << " doubles -> to FITS file" << dtfile << endl;
218 xntms_=new double[nchan*nband+2];
219}
220
221/* --Methode-- */
222void BRMeanSpecCalculator::SetNumberOfBands(int numberOfBands, int ibandfirst, int ibandlast)
223{
224 if (numberOfBands < 1) numberOfBands = 1;
225 numberOfBands_ = numberOfBands;
226 if (ibandfirst < 0 )ibandfirst = 0;
227 if (ibandlast >= numberOfBands_ ) ibandlast = numberOfBands_-1;
228 ibandfirst_=ibandfirst; ibandlast_=ibandlast;
229
230 cout << "BRMeanSpecCalculator::SetNumberOfBands (END) : "
231 << numberOfBands_ << " "
232 << ibandfirst_ << " "
233 << ibandlast_ << endl;
234
235}
236
237/* --Methode-- */
238void BRMeanSpecCalculator::ReadGainFitsFile(string filename, bool fgapp)
239{
240 cout << " BRMeanSpecCalculator::ReadGainFitsFile() - reading file " << filename;
241 FitsInOutFile fis(filename, FitsInOutFile::Fits_RO);
242 fis >> sgain_;
243 fg_apply_gains_=fgapp;
244 cout << " MeanGain=" << sgain_.Sum()/sgain_.Size() << " ApplyGains="
245 << ((fg_apply_gains_)?"true":"false") << endl;
246 if( (spec_window_.SizeX()!= sgain_.NCols()) || (spec_window_.SizeY()!= sgain_.NRows()) ){
247 cout << " ReadGainFitsFile: BAD Gain Matrix sizes " << endl;
248 sgain_.Show();
249 spec_window_.Show();
250 throw ParmError("ReadGainFitsFile: BAD Gain Matrix sizes");
251 }
252}
253
254//JEC
255//static inline r_4 Zmod2(complex<r_4> z)
256//{ return (z.real()*z.real()+z.imag()*z.imag()); }
257
258
259
260
261/* --Methode-- */
262int BRMeanSpecCalculator::Process()
263{
264 // Cette methode remplit le tableau spec_window_ avec les spectres (renormalise avec
265 // les gains si demande) et appelle la methode du traitement de la fenetre temporelle
266 // des spectres le cas echeant ProcSpecWin()
267
268
269 int_8 nbpaqdec = (int_8)totnbpaq_-(int_8)GetSpecWinOverlapSize();
270 if ((nbpaqdec>0)&&(nbpaqdec%GetSpectraWindowSize()==0)) {
271 paqnum_w_end=totnbpaq_-GetSpecWinExtensionSize();
272 ProcSpecWin(paqnum_w_start, paqnum_w_end);
273 paqnum_w_start=totnbpaq_-GetSpecWinExtensionSize();
274 }
275
276 if (fgdatafft_) { // Donnees firmware FFT
277 for(sa_size_t i=0; i<spec_window_.SizeY(); i++) {
278 TwoByteComplex* zp=NULL;
279 if (fgsinglechannel_) {
280 zp=vpaq_[i].Data1C();
281 }
282 else {
283 zp=vpaq_[i/2].Data1C();
284 if (i%2==1) zp=vpaq_[i/2].Data2C();
285 }
286 sa_size_t kz=PaqNumToArrayIndex(totnbpaq_);
287 for(sa_size_t f=0; f<spec_window_.SizeX(); f++)
288 spec_window_(f,i,kz) = zp[f].module2F();
289 }
290 }
291 else { // Donnees RAW qui ont du etre processe par BRFFTCalculator
292 for(sa_size_t i=0; i<spec_window_.SizeY(); i++) {
293 complex<ODT>* zp=NULL;
294 if (fgsinglechannel_) {
295 zp=reinterpret_cast< complex<ODT>* > (vprocpaq_[i]);
296 }
297 else {
298 zp=reinterpret_cast< complex<ODT>* > (vprocpaq_[i/2]);
299 if (i%2==1) zp= reinterpret_cast< complex<ODT>* >(vprocpaq_[i/2]+memgr_.ProcPaqSize()/2) ;
300 }
301 sa_size_t kz=PaqNumToArrayIndex(totnbpaq_);
302 for(sa_size_t f=0; f<spec_window_.SizeX(); f++)
303 spec_window_(f,i,kz) = Zmod2(zp[f]);
304 }
305 }
306 if (fg_apply_gains_) { // Application des gains, si demande
307 sa_size_t kz=PaqNumToArrayIndex(totnbpaq_);
308 for(sa_size_t i=0; i<spec_window_.SizeY(); i++) {
309 (spec_window_(Range::all(), Range(i), Range(kz)).CompactAllDimensions()).Div(sgain_.Row(i).CompactAllDimensions());
310 }
311 }
312
313 // Pour calcul MeanTimeTag , MeanFrameCounter
314 moyfc_+=curfc_[0];
315 moytt_+=((double)vpaq_[0].TimeTag()/1.25e8);
316 nbpmoyttfc_++;
317 totnbpaq_++;
318 return 0;
319}
320
321
322/* --Methode-- */
323void BRMeanSpecCalculator::ProcSpecWin(uint_8 numpaqstart, uint_8 numpaqend)
324{
325
326 if (prtlev_>0) {
327 uint_8 modulo = prtmodulo_/GetSpectraWindowSize();
328 if (modulo<1) modulo=1;
329 if (nbtot_specwin_%modulo==0) {
330 cout << " BRMeanSpecCalculator::ProcSpecWin() num_win=" << nbtot_specwin_ << " numpaqstart=" << numpaqstart
331 << " numpaqend=" << numpaqend << endl;
332 cout << " ... ObsTime=" << getObsTime() << " TimeTag=" << getCurTimeTagSeconds() << " s. FrameCounter="
333 << getCurFrameCounter() << endl;
334 }
335 }
336 // On appelle la routine de nettoyage qui doit flagger les mauvais paquets
337 FlagBadPackets(numpaqstart, numpaqend);
338
339 // Boucle sur les numeros de paquets de la fenetre en temps
340 for (uint_8 jp=numpaqstart; jp<numpaqend; jp++) {
341 // On sauvegarde les spectres moyennes si necessaire
342 if ((nbpaq4mean_[0]>0)&&(nbpaq4mean_[0]%nmean_ == 0)) SaveMeanSpectra();
343 // On peut aussi acceder aux spectres et flags pour (jpmin -),(jpmax+) GetSpecWinExtensionSize()
344 sa_size_t kz=PaqNumToArrayIndex(jp);
345 // Boucle sur les numeros de voie (canaux)
346 for(sa_size_t i=0; i<spec_window_.SizeY(); i++) {
347 if ( clnflg_(i,kz) != 0) continue;
348 TVector< r_4 > spec = mspecmtx_.Row(i);
349 TVector< r_4 > sspec = sigspecmtx_.Row(i);
350 // Calcul de spectres moyennes et variance
351 for(sa_size_t f=1; f<spec.Size(); f++) { // boucle sur les frequences
352 spec(f) += spec_window_(f,i,kz);
353 sspec(f) += spec_window_(f,i,kz)*spec_window_(f,i,kz);
354 }
355 nbpaq4mean_[i]++; // compteur de paquets OK pour la moyenne
356 }
357 }
358 if (nbtot_specwin_<nmaxfiles_specw_) SaveSpectraWindow();
359 nbtot_specwin_++;
360 return;
361}
362
363/* --Methode-- */
364void BRMeanSpecCalculator::FlagBadPackets(uint_8 numpaqstart, uint_8 numpaqend)
365{
366 // Boucle sur les numeros de paquets de la fenetre en temps
367 for (uint_8 jp=numpaqstart; jp<numpaqend; jp++) {
368 // On peut aussi acceder aux spectres et flags pour (jpmin -),(jpmax+) GetSpecWinExtensionSize()
369 sa_size_t kz=PaqNumToArrayIndex(jp);
370 // Boucle sur les numeros de voie (canaux)
371 for(sa_size_t i=0; i<spec_window_.SizeY(); i++) {
372
373 clnflg_(i,kz)=0;
374
375 double mean, sigma;
376 ////////BUG sa_size_t kz=PaqNumToArrayIndex(totnbpaq_);
377 double variance=0.;
378 variance=spec_window_(Range(1,Range::lastIndex()), Range(i), Range(kz)).Sum();
379 if(xnt_)xnt_[i]=variance;
380 if(numberOfBands_>0){
381 //Compute nomalized variance in bands freq.
382 sa_size_t fMin;
383 sa_size_t fMax;
384 int bandW = spec_window_.SizeX()/numberOfBands_;
385 vector<double> varNomBinned(numberOfBands_);
386 for (sa_size_t j=ibandfirst_; j<=ibandlast_; j++){
387 fMin = j*bandW;
388 fMax =fMin+bandW-1;
389 varNomBinned[j]=spec_window_(Range(fMin,fMax), Range(i), Range(kz)).Sum();
390 varNomBinned[j]/=(r_4)bandW;
391 if(xnt_)xnt_[spec_window_.SizeY()+i*numberOfBands_+j] = varNomBinned[j];
392
393 // cout << "(jec) var["<<j<<"] =" << varNomBinned[j]
394 // << " min " << varmin_
395 // << " max " << varmax_ << endl;
396 if(varNomBinned[j]<varmin_) { clnflg_(i,kz)=10+j; nbadpaq_[i]++; break;}
397 else if(varNomBinned[j]>varmax_) { clnflg_(i,kz)=100+j; nbadpaq_[i]++; break;}
398 }
399 //cout << "clnflg_("<<i<<","<<kz<<"): " << clnflg_(i,kz) << endl;
400 }//if bands
401 }//loop on channels
402
403 if (dtp_ && xnt_) dtp_->AddRow(xnt_);
404 }
405 return;
406}
407
408/* --Methode-- */
409void BRMeanSpecCalculator::SaveMeanSpectra()
410{
411 for(sa_size_t ir=0; ir<mspecmtx_.NRows(); ir++) {
412 char buff[32];
413 sprintf(buff,"NPAQSUM_%d",(int)ir);
414 mspecmtx_.Info()["NPAQSUM"] = nbpaq4mean_[0];
415 mspecmtx_.Info()[buff] = nbpaq4mean_[ir];
416 sigspecmtx_.Info()["NPAQSUM"] = nbpaq4mean_[0];
417 sigspecmtx_.Info()[buff] = nbpaq4mean_[ir];
418 if (nbpaq4mean_[ir] > 0) {
419 mspecmtx_.Row(ir) /= (r_4)nbpaq4mean_[ir];
420 sigspecmtx_.Row(ir) /= (r_4)nbpaq4mean_[ir];
421 sigspecmtx_.Row(ir) -= (mspecmtx_.Row(ir) && mspecmtx_.Row(ir)); // Mean(X^2) - [ Mean(X) ]^2
422 }
423 }
424 mspecmtx_.Info()["ENDTIME"]=getObsTime();
425 sigspecmtx_.Info()["ENDTIME"]=getObsTime();
426 if (nbpmoyttfc_>1) {
427 moyfc_/=(double)nbpmoyttfc_;
428 moytt_/=(double)nbpmoyttfc_;
429 }
430 string ikey,ikdesc;
431 ikey="MeanFC"; ikdesc="Mean FrameCounter";
432 mspecmtx_.Info().SetI(ikey,moyfc_);
433 mspecmtx_.Info().SetComment(ikey,ikdesc);
434 sigspecmtx_.Info().SetI(ikey,moyfc_);
435 sigspecmtx_.Info().SetComment(ikey,ikdesc);
436
437 ikey="MeanTT"; ikdesc="Mean TimeTag";
438 mspecmtx_.Info().SetD(ikey,moytt_);
439 mspecmtx_.Info().SetComment(ikey,ikdesc);
440 sigspecmtx_.Info().SetD(ikey,moytt_);
441 sigspecmtx_.Info().SetComment(ikey,ikdesc);
442
443 char nfile[64];
444 string flnm;
445 {
446 sprintf(nfile,"mspecmtx%d.fits",numfile_);
447 flnm="!"+outpath_+nfile;
448 FitsInOutFile fos(flnm,FitsInOutFile::Fits_Create);
449 fos << mspecmtx_;
450 }
451 {
452 sprintf(nfile,"sigspecmtx%d.fits",numfile_);
453 flnm="!"+outpath_+nfile;
454 FitsInOutFile fos(flnm,FitsInOutFile::Fits_Create);
455 fos << sigspecmtx_;
456 }
457
458 cout << numfile_ << "-BRMeanSpecCalculator::SaveMeanSpectra() NPaqProc="
459 << totnbpaq_ << " -> Mean/Sig spectra Matrix in " << flnm << " /sigspec...ppf" << endl;
460 numfile_++;
461
462 if (dtpms_!=NULL) FillPwrTmDTable(mspecmtx_);
463
464 for(size_t i=0; i<nbpaq4mean_.size(); i++) nbpaq4mean_[i]=0;
465 mspecmtx_ = (r_4)(0.);
466 sigspecmtx_ = (r_4)(0.);
467 moyfc_=moytt_=0.;
468 nbpmoyttfc_=0;
469 return;
470}
471
472/* --Methode-- */
473void BRMeanSpecCalculator::SaveSpectraWindow()
474{
475 char nfile[64];
476 string flnm;
477 sprintf(nfile,"specwin%d.fits",nbtot_specwin_);
478 flnm="!"+outpath_+nfile;
479 FitsInOutFile fos(flnm,FitsInOutFile::Fits_Create);
480 fos << spec_window_;
481 cout << " SaveSpectraWindow() " << nbtot_specwin_ << "- file " << nfile << " created " << endl;
482}
483
484/* --Methode-- */
485void BRMeanSpecCalculator::FillPwrTmDTable(TMatrix< r_4 >& specmtx)
486{
487 if (dtpms_==NULL) return;
488 xntms_[0]=moytt_;
489 xntms_[1]=getObsTimeSeconds();
490 size_t knt=2;
491 sa_size_t wband=(lastfreqms_-firstfreqms_)/nbandms_;
492 for(sa_size_t ir=0; ir<specmtx.NRows(); ir++) {
493 for(sa_size_t jf=0; jf<specmtx.NCols()-wband; jf+=wband) {
494 TVector<r_4> srow=specmtx.Row(ir);
495 xntms_[knt]=srow.SubVector(Range(jf,jf+wband-1)).Sum();
496 knt++;
497 }
498 }
499 dtpms_->AddRow(xntms_);
500}
501
502//---------------------------------------------------------------------
503// Classe de thread de calcul de FFT sur donnees RAW
504//---------------------------------------------------------------------
505/*!
506 \class BRFFTCalculator
507 \ingroup TAcq
508
509 \brief Computes FFT on data from RAW ADC board firmware.
510*/
511/* --Methode-- */
512BRFFTCalculator::BRFFTCalculator(RAcqMemZoneMgr& memgr, bool fgsinglechannel)
513 : BRBaseProcessor(memgr), fgsinglechannel_(fgsinglechannel), totnbfftpaq_(0)
514{
515 uint_4 nb_octets_entrop = 0; //this value is valid for Dec. 2010 data at Nancay
516 const char* venvp = NULL;
517 venvp=getenv("BRANA_NBYTECUT");
518 if (venvp!=NULL){
519 nb_octets_entrop = atoi(venvp);
520 cout << "BRFFTCalculator : BRANA_NBYTECUT : " << nb_octets_entrop << endl;
521 }
522
523 BRPaquet paq(memgr_.PaqSize()-nb_octets_entrop);
524 //JEC END
525 // BRPaquet paq(memgr_.PaqSize());
526 setNameId("FFTCalc",2);
527 ffts_.SetInDataSize((fgsinglechannel_)?paq.DataSize():paq.DataSize()/2);
528}
529
530/* --Methode-- */
531BRFFTCalculator::~BRFFTCalculator()
532{
533}
534
535
536/* --Methode-- */
537int BRFFTCalculator::Process()
538{
539 for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
540 ffts_.DoFFT( reinterpret_cast<IDT *>(vpaq_[fib].Data1() ),
541 reinterpret_cast< complex<ODT>* > (vprocpaq_[fib]) );
542 totnbfftpaq_++;
543 if ( fgsinglechannel_ ) continue;
544 ffts_.DoFFT( reinterpret_cast<IDT *>(vpaq_[fib].Data2() ),
545 reinterpret_cast< complex<ODT>* > (vprocpaq_[fib]+memgr_.ProcPaqSize()/2) );
546 totnbfftpaq_++;
547 }
548 return 0;
549}
550
551
552//-------------------------------------------------------------------------
553// Classe WBRFFT : Calcul de TF sur donnees brutes (firmware RAW)
554//-------------------------------------------------------------------------
555/*!
556 \class WBRFFT
557 \ingroup TAcq
558
559 \brief Helper class for FFT computation used by BRFFTCalculator class.
560
561 This class uses FFTW library.
562*/
563
564ZMutex* WBRFFT::mtx_fftwp_ = NULL;
565
566/* --Methode-- */
567WBRFFT::WBRFFT(uint_4 sz)
568 : sz_(sz)
569{
570 if (mtx_fftwp_ == NULL) mtx_fftwp_ = new ZMutex;
571 if (sz>0) SetInDataSize(sz);
572}
573
574/* --Methode-- */
575WBRFFT::~WBRFFT()
576{
577}
578
579/* --Methode-- */
580void WBRFFT::SetInDataSize(uint_4 sz)
581{
582 sz_ = sz;
583 if (sz_<1) return;
584 inp.SetSize(sz);
585 outfc.SetSize(sz/2+1);
586 mtx_fftwp_->lock();
587 myplan_ = fftwf_plan_dft_r2c_1d(inp.Size(), inp.Data(),
588 (fftwf_complex*)outfc.Data(), FFTW_ESTIMATE);
589 mtx_fftwp_->unlock();
590}
591
592/* --Methode-- */
593void WBRFFT::DoFFT( IDT *indata, complex<ODT> * ofc)
594{
595 if (sz_<1) return;
596 for(uint_4 k=0; k<inp.Size(); k++) inp(k)=(ODT)indata[k];
597 fftwf_execute(myplan_);
598 for(uint_4 k=0; k<outfc.Size(); k++) ofc[k]=outfc(k)/(ODT)sz_; // on renormalise les coeff FFT ( / sz )
599 return;
600}
601
602/* --Methode-- */
603void WBRFFT::PrintData(IDT *indata, complex<ODT> * ofc, uint_4 sz)
604{
605 cout << " --- WBRFFT::PrintData() size=" << sz << endl;
606 for(uint_4 k=0; k<sz; k+=8) {
607 IDT* in = indata+k;
608 cout << " Indata[" << k << "..." << k+8 << "]= ";
609 for(uint_4 i=0; i<8; i++) cout << (IIDT)in[i] << " ";
610 cout << endl;
611 }
612 cout << endl;
613 for(uint_4 k=0; k<sz/2; k+=4) {
614 complex< ODT>* out = ofc+k;
615 cout << " OutFC[" << k << "..." << k+4 << "]= ";
616 for(uint_4 i=0; i<4; i++) cout << out[i] << " ";
617 cout << endl;
618 }
619 return;
620
621}
622
623
624//---------------------------------------------------------------
625// Classe thread de traitement donnees ADC avec 2 voies par frame
626// !!!! OBSOLETE !!!!
627//---------------------------------------------------------------
628/*!
629 \class BRProcA2C
630 \ingroup TAcq
631
632 \brief Deprecated processing class.
633*/
634
635// Mutex pour eviter le plantage du a FFTW qui ne semble pas thread-safe
636static ZMutex* pmutfftw=NULL;
637
638/* --Methode-- */
639BRProcA2C::BRProcA2C(RAcqMemZoneMgr& mem, string& path, bool fgraw, uint_4 nmean,
640 uint_4 nmax, bool fghist, uint_4 nfsmap, bool fgnotrl, int card)
641 : memgr(mem)
642{
643 fgraw_ = fgraw;
644 nmax_ = nmax;
645 nmean_ = nmean;
646 if (fgraw_) cout << " BRProcA2C::BRProcA2C() - constructeur RAW data - NMean=" << nmean_ << endl;
647 else cout << " BRProcA2C::BRProcA2C() - constructeur FFT data - NMean=" << nmean_ << endl;
648 nfsmap_ = nfsmap;
649 stop_ = false;
650 path_ = path;
651 fgnotrl_ = fgnotrl;
652 fghist_ = fghist;
653 card_ = card;
654 if (pmutfftw==NULL) pmutfftw=new ZMutex;
655}
656
657/* --Methode-- */
658void BRProcA2C::Stop()
659{
660 stop_=true;
661 // cout <<" BRProcA2C::Stop ... > STOP " << endl;
662}
663
664
665
666
667static inline string card2name_(int card)
668{
669 if (card==2) return " (Chan3,4) ";
670 else return " (Chan1,2) ";
671}
672/* --Methode-- */
673void BRProcA2C::run()
674{
675 setRC(1);
676 try {
677 Timer tm("BRProcA2C", false);
678 TimeStamp ts;
679 BRPaqChecker pcheck(!fgnotrl_); // Verification/comptage des paquets
680
681 size_t totnbytesout = 0;
682 size_t totnbytesproc = 0;
683
684 cout << " BRProcA2C::run() - Starting " << ts << " NMaxMemZones=" << nmax_
685 << " NMean=" << nmean_ << card2name_(card_) << endl;
686 cout << " BRProcA2C::run()... - Output Data Path: " << path_ << endl;
687 char fname[512];
688// sprintf(fname,"%s/proc.log",path_.c_str());
689// ofstream filog(fname);
690// filog << " BRProcA2C::run() - starting log file " << ts << endl;
691// filog << " ... NMaxMemZones=" << nmax_ << " NMean=" << nmean_ << " Step=" << step_ << endl;
692
693/*----DELETE NTuple
694 const char* nnames[8] = {"fcs","tts","s1","s2","s12","s12re","s12im","s12phi"};
695 NTuple nt(8, nnames);
696 double xnt[10];
697 uint_4 nmnt = 0;
698 double ms1,ms2,ms12,ms12re,ms12im,ms12phi;
699----*/
700// Time sample (raw data) /FFT coeff histograms
701 Histo* ph1=NULL;
702 Histo* ph2=NULL;
703 if (fghist_) {
704 if (fgraw_) {
705 ph1 = new Histo(-0.5, 255.5, 256);
706 ph2 = new Histo(-0.5, 255.5, 256);
707 }
708 else {
709 ph1 = new Histo(-128.5, 128.5, 257);
710 ph2 = new Histo(-128.5, 128.5, 257);
711 }
712 }
713
714// Initialisation pour calcul FFT
715 TVector< complex<r_4> > cfour1; // composant TF
716 uint_4 paqsz = memgr.PaqSize();
717 uint_4 procpaqsz = memgr.ProcPaqSize();
718
719
720 BRPaquet pq(NULL, NULL, paqsz);
721 TVector<r_4> vx(pq.DataSize()/2);
722 int szfour = pq.DataSize()/2/2+1;
723 cfour1.SetSize(szfour);
724/*
725 vx = (r_4)(0.);
726 FFTPackServer ffts;
727 ffts.FFTForward(vx, cfour1);
728 szfour = cfour1.Size();
729*/
730
731 bool fgtimfreq = false; // true->cartes temps<>frequences
732 if (nfsmap_>0) fgtimfreq=true;
733
734 TVector< complex<r_4> > cfour2(cfour1.Size());
735
736 TVector<r_4> spectreV1(cfour1.Size());
737 TVector<r_4> spectreV2(cfour1.Size());
738 TVector<r_4> moyspecV1(cfour1.Size()); // Moyenne des Spectres
739 TVector<r_4> moyspecV2(cfour1.Size());
740 TVector<r_4> sigspecV1(cfour1.Size()); // Sigma des Spectres
741 TVector<r_4> sigspecV2(cfour1.Size());
742 TVector< complex<r_4> > visiV12( cfour1.Size() );
743
744 TMatrix<r_4> timfreqV1, timfreqV2; // Cartes temps<>frequences
745 if (fgtimfreq) {
746 timfreqV1.SetSize(nmean_, spectreV1.Size()/nfsmap_);
747 timfreqV2.SetSize(nmean_, spectreV2.Size()/nfsmap_);
748 }
749 cout << " *DBG*BRProcA2C PaqSz=" << paqsz << " ProcPaqSize=" << procpaqsz
750 << " procpaqsz/2=" << procpaqsz/2 << " cfour1.Size()=" << cfour1.Size()
751 << " *8=" << cfour1.Size()*8 << endl;
752
753 pmutfftw->lock();
754 fftwf_plan plan1 = fftwf_plan_dft_r2c_1d(vx.Size(), vx.Data(),
755 (fftwf_complex*)cfour1.Data(), FFTW_ESTIMATE);
756 fftwf_plan plan2 = fftwf_plan_dft_r2c_1d(vx.Size(), vx.Data(),
757 (fftwf_complex*)cfour2.Data(), FFTW_ESTIMATE);
758 pmutfftw->unlock();
759
760 uint_4 ifile = 0;
761 uint_4 nzm = 0; // Nb de paquets moyennes pour le calcul de chaque spectre
762 uint_4 nmoyspec = 0; // Nb de spectres moyennes
763
764 uint_4 curfc=0;
765 uint_8 curtt=0;
766 uint_8 firsttt=0;
767 bool fgfirst=true;
768 double moysig[2]={0.,0.};
769 double sigsig[2]={0.,0.};
770 uint_8 nbsig[2]={0,0};
771
772 for (uint_4 kmz=0; kmz<nmax_; kmz++) {
773 if (stop_) break;
774 int mid = memgr.FindMemZoneId(MemZA_ProcA);
775 Byte* buff = memgr.GetMemZone(mid);
776 if (buff == NULL) {
777 cout << " BRProcA2C::run()/ERROR memgr.GetMemZone(" << mid << ") -> NULL" << endl;
778 break;
779 }
780 Byte* procbuff = memgr.GetProcMemZone(mid);
781 if (procbuff == NULL) {
782 cout << " BRProcA2C::run()/ERROR memgr.GetProcMemZone(" << mid << ") -> NULL" << endl;
783 break;
784 }
785//---- DELETE nmnt=0; ms1=ms2=ms12=ms12re=ms12im=ms12phi=0.;
786 for(uint_4 i=0; i<memgr.NbPaquets(); i++) {
787 BRPaquet paq(NULL, buff+i*paqsz, paqsz);
788 if (!pcheck.Check(paq)) continue; // on ne traite que les paquets OK
789 if (fgfirst) { firsttt=paq.TimeTag(); fgfirst=false; }
790 curfc=paq.FrameCounter();
791 curtt=paq.TimeTag()-firsttt;
792// Traitement voie 1
793 if (fghist_) {
794 for(sa_size_t j=0; j<vx.Size(); j++) {
795 r_4 vts=(fgraw_)?((r_4)(*(paq.Data1()+j))):((r_4)(*(paq.Data1S()+j)));
796 ph1->Add((r_8)vts);
797 moysig[0] += (double)vts;
798 sigsig[0] += ((double)vts)*((double)vts);
799 nbsig[0]++;
800 }
801 for(sa_size_t j=0; j<vx.Size(); j++) {
802 r_4 vts=(fgraw_)?((r_4)(*(paq.Data2()+j))):((r_4)(*(paq.Data2S()+j)));
803 ph2->Add((r_8)vts);
804 moysig[1] += (double)vts;
805 sigsig[1] += ((double)vts)*((double)vts);
806 nbsig[1]++;
807 }
808 }
809 if (fgraw_) {
810 for(sa_size_t j=0; j<vx.Size(); j++)
811 vx(j) = (r_4)(*(paq.Data1()+j))-127.5;
812 // fftwf_complex* coeff1 = (fftwf_complex*)(procbuff+i*procpaqsz);
813 fftwf_execute(plan1);
814 // Traitement voie 2
815 for(sa_size_t j=0; j<vx.Size(); j++)
816 vx(j) = (r_4)(*(paq.Data2()+j))-127.5;
817 fftwf_execute(plan2);
818 }
819 else {
820 for(sa_size_t j=1; j<cfour1.Size()-1; j++) {
821 cfour1(j) = complex<r_4>((r_4)paq.Data1C()[j].realB(), (r_4)paq.Data1C()[j].imagB());
822 cfour2(j) = complex<r_4>((r_4)paq.Data2C()[j].realB(), (r_4)paq.Data2C()[j].imagB());
823 }
824 cfour1(0) = complex<r_4>((r_4)paq.Data1C()[0].realB(), (r_4)0.);
825 cfour1(cfour1.Size()-1) = complex<r_4>((r_4)paq.Data1C()[0].imagB(), (r_4)0.);
826 cfour2(0) = complex<r_4>((r_4)paq.Data2C()[0].realB(), (r_4)0.);
827 cfour2(cfour2.Size()-1) = complex<r_4>((r_4)paq.Data2C()[0].imagB(), (r_4)0.);
828 }
829 for(sa_size_t j=0; j<spectreV1.Size(); j++)
830 spectreV1(j) += BRMeanSpecCalculator::Zmod2(cfour1(j));
831 memcpy(procbuff+i*procpaqsz, cfour1.Data(), sizeof(complex<r_4>)*cfour1.Size());
832 if (fgtimfreq) { // Remplissage tableau temps-frequence
833 for(sa_size_t c=1; c<timfreqV1.NCols(); c++) {
834 for(sa_size_t j=c*nfsmap_; j<(c+1)*nfsmap_; j++)
835 timfreqV1(nzm, c) += BRMeanSpecCalculator::Zmod2(cfour1(j));
836 }
837 }
838 for(sa_size_t j=0; j<spectreV2.Size(); j++)
839 spectreV2(j) += BRMeanSpecCalculator::Zmod2(cfour2(j)); // BRMeanSpecCalculator::Zmod2(zp2[j]);
840 memcpy(procbuff+i*procpaqsz+procpaqsz/2, cfour2.Data(), sizeof(complex<r_4>)*cfour2.Size());
841 if (fgtimfreq) { // Remplissage tableau temps-frequence
842 for(sa_size_t c=1; c<timfreqV2.NCols(); c++) {
843 for(sa_size_t j=c*nfsmap_; j<(c+1)*nfsmap_; j++)
844 timfreqV2(nzm,c) += BRMeanSpecCalculator::Zmod2(cfour2(j));
845 }
846 }
847
848// Calcul correlation (visibilite V1 * V2)
849 for(sa_size_t j=0; j<visiV12.Size(); j++)
850 visiV12(j)+=cfour1(j)*conj(cfour2(j));
851// for(sa_size_t j=0; j<visiV12.Size(); j++) visiV12(j)+=zp1[j]*zp2[j];
852 if (nzm==0) {
853 spectreV1.Info()["StartFC"] = curfc;
854 spectreV2.Info()["StartFC"] = curfc;
855 visiV12.Info()["StartFC"] = curfc;
856 spectreV1.Info()["StartTT"] = curtt;
857 spectreV2.Info()["StartTT"] = curtt;
858 visiV12.Info()["StartTT"] = curtt;
859 }
860 nzm++;
861/*----DELETE
862 if (nmnt==0) { xnt[0]=paq.FrameCounter(); xnt[1]=paq.TimeTag(); }
863 for(sa_size_t j=2700; j<2800; j++) {
864 ms1 += BRMeanSpecCalculator::Zmod2(cfour1(j)); ms2 += BRMeanSpecCalculator::Zmod2(cfour2(j));
865 complex<r_4> zvis = cfour1(j)*conj(cfour2(j));
866 ms12 += BRMeanSpecCalculator::Zmod2(zvis); ms12re += zvis.real(); ms12im += zvis.imag();
867 ms12phi+= atan2(zvis.imag(),zvis.real());
868 }
869 nmnt++;
870----*/
871 totnbytesproc += paq.DataSize();
872 totnbytesout += (2*sizeof(complex<r_4>)*cfour1.Size());
873
874 } // Fin de boucle sur les paquets d'une zone
875
876/*---- DELETE
877 if (nmnt>0) {
878 double fnorm = (double)nmnt*(2800-2700);
879 xnt[2] = ms1 /= fnorm;
880 xnt[3] = ms2 /= fnorm;
881 xnt[4] = ms12 /= fnorm;
882 xnt[5] = ms12re /= fnorm;
883 xnt[6] = ms12im /= fnorm;
884 xnt[7] = ms12phi /= fnorm;
885 nt.Fill(xnt);
886 }
887----*/
888 if ((nzm >= nmean_) || ((kmz==(nmax_-1))&&(nzm>1))) {
889 spectreV1 /= (r_4)(nzm);
890 spectreV2 /= (r_4)(nzm);
891
892 // pour le calcul des moyennes et sigmas de ces spectres
893 moyspecV1 += spectreV1;
894 moyspecV2 += spectreV2;
895 sigspecV1 += (spectreV1 && spectreV1);
896 sigspecV2 += (spectreV2 && spectreV2);
897 nmoyspec++;
898
899 visiV12 /= complex<r_4>((r_4)nzm, 0.);
900
901 spectreV1.Info()["NPaqMoy"] = nzm;
902 spectreV2.Info()["NPaqMoy"] = nzm;
903 visiV12.Info()["NPaqMoy"] = nzm;
904 spectreV1.Info()["EndFC"] = curfc;
905 spectreV2.Info()["EndFC"] = curfc;
906 visiV12.Info()["EndFC"] = curfc;
907 spectreV1.Info()["EndTT"] = curtt;
908 spectreV2.Info()["EndTT"] = curtt;
909 visiV12.Info()["EndTT"] = curtt;
910 {
911 sprintf(fname,"%s_%d.ppf",path_.c_str(),(int)ifile);
912 POutPersist po(fname);
913 string tag1="specV1";
914 string tag2="specV2";
915 string tag12="visiV12";
916 string tagh1="tshV1";
917 string tagh2="tshV2";
918 string tagtf1="timfreqV1";
919 string tagtf2="timfreqV2";
920 if (card_==2) {
921 tag1 = "specV3";
922 tag2 = "specV4";
923 tagh1 = "tshV1";
924 tagh2 = "tshV2";
925 tag12="visiV34";
926 tagtf1="timfreqV3";
927 tagtf2="timfreqV4";
928 }
929 po << PPFNameTag(tag1) << spectreV1;
930 po << PPFNameTag(tag2) << spectreV2;
931 po << PPFNameTag(tag12) << visiV12;
932 if (fghist_) {
933 po << PPFNameTag(tagh1) << (*ph1);
934 po << PPFNameTag(tagh2) << (*ph2);
935
936 double sspvmax[3] = {0.,0.,0.};
937 int_4 sspvmaxidx[3] = {-1,-1,-1};
938 for(int jji=1;jji<visiV12.Size()-1;jji++) {
939 r_4 zmv2 = BRMeanSpecCalculator::Zmod2(visiV12(jji));
940 if (zmv2>sspvmax[2]) { sspvmax[2]=zmv2; sspvmaxidx[2]=jji; }
941 }
942 TVector<r_4>& sspv = spectreV1;
943 for(int ic=0; ic<2; ic++) {
944 if (ic==1) sspv = spectreV2;
945 for(int jji=1;jji<sspv.Size()-1;jji++)
946 if (sspv(jji)>sspvmax[ic]) { sspvmax[ic]=sspv(jji); sspvmaxidx[ic]=jji; }
947 if (nbsig[ic] < 1) { moysig[ic]=sigsig[ic]=-1.; }
948 else {
949 moysig[ic] /= (double)nbsig[ic];
950 sigsig[ic] /= (double)nbsig[ic];
951 sigsig[ic] -= (moysig[ic]*moysig[ic]);
952 sigsig[ic] = sqrt(sigsig[ic]);
953 cout << "===Voie " << ic << " Moy=" << moysig[ic] << " Sig=" << sigsig[ic]
954 << " MaxSpec Amp= " << sqrt(sspvmax[ic])/double(pq.DataSize()/2/2)
955 << " Pos=" << sspvmaxidx[ic] << " (NPts=" << nbsig[ic] << ")" << endl;
956 }
957 }
958 cout << "=== Voie1x2 MaxSpec Amp= " << sqrt(sqrt(sspvmax[2])/double(pq.DataSize()/2/2))
959 << " Pos=" << sspvmaxidx[2] << endl;
960 } // fin if (fghist_)
961
962 if (fgtimfreq) {
963 timfreqV1 /= (r_4)nzm;
964 timfreqV2 /= (r_4)nzm;
965 po << PPFNameTag(tagtf1) << timfreqV1;
966 po << PPFNameTag(tagtf2) << timfreqV2;
967 }
968 }
969 spectreV1 = (r_4)(0.);
970 spectreV2 = (r_4)(0.);
971 visiV12 = complex<r_4>(0., 0.);
972 if (fghist_) {
973 ph1->Zero();
974 ph2->Zero();
975 moysig[0]=moysig[1]=0.;
976 sigsig[0]=sigsig[1]=0.;
977 nbsig[0]=nbsig[1]=0;
978 }
979 if (fgtimfreq) {
980 timfreqV1 = (r_4)(0.);
981 timfreqV2 = (r_4)(0.);
982 }
983 nzm = 0; ifile++;
984// ts.SetNow();
985// filog << ts << " : proc file " << fname << endl;
986 cout << " BRProcA2C::run() created file " << fname << card2name_(card_) << endl;
987 }
988
989 memgr.FreeMemZone(mid, MemZS_ProcA);
990 } // Fin de boucle sur les zones a traiter
991 cout << " ------------ BRProcA2C::run() END " << card2name_(card_)
992 << " ------------ " << endl;
993/*---- DELETE
994 {
995 nt.Info()["FirstTT"]=firsttt;
996 cout << nt;
997 sprintf(fname,"%s_nt.ppf",path_.c_str());
998 POutPersist po(fname);
999 po << PPFNameTag("ntv12") << nt;
1000 cout << " BRProcA2C::run() created NTuple file " << fname << card2name_(card_) << endl;
1001 }
1002---- */
1003 if (nmoyspec>0) { // Calcul des moyennes et sigmas des spectres
1004 r_4 fnms = nmoyspec;
1005 moyspecV1 /= fnms;
1006 moyspecV2 /= fnms;
1007 sigspecV1 /= fnms;
1008 sigspecV2 /= fnms;
1009 sigspecV1 -= (moyspecV1 && moyspecV1);
1010 sigspecV2 -= (moyspecV2 && moyspecV2);
1011 sigspecV1 = Sqrt(sigspecV1);
1012 sigspecV2 = Sqrt(sigspecV2);
1013 TVector<r_4> rsbV1, rsbV2; // Rapport signal/bruit
1014 moyspecV1.DivElt(sigspecV1, rsbV1, false, true);
1015 moyspecV2.DivElt(sigspecV2, rsbV2, false, true);
1016 sprintf(fname,"%s_ms.ppf",path_.c_str());
1017 POutPersist po(fname);
1018 po << PPFNameTag("moyspecV1") << moyspecV1;
1019 po << PPFNameTag("moyspecV2") << moyspecV2;
1020 po << PPFNameTag("sigspecV1") << sigspecV1;
1021 po << PPFNameTag("sigspecV2") << sigspecV2;
1022 po << PPFNameTag("rsbV1") << rsbV1;
1023 po << PPFNameTag("rsbV2") << rsbV2;
1024 cout << " BRProcA2C::run() created moysigspec file " << fname << card2name_(card_) << endl;
1025 }
1026
1027 if (fghist_) {
1028 delete ph1;
1029 delete ph2;
1030 }
1031 ts.SetNow();
1032 tm.SplitQ();
1033 cout << " TotalProc= " << totnbytesproc/(1024*1024) << " MBytes, rate= "
1034 << (double)(totnbytesproc)/1024./tm.PartialElapsedTimems() << " MB/s"
1035 << " ProcDataOut=" << totnbytesout/(1024*1024) << " MB" << endl;
1036 cout << pcheck;
1037 cout << " BRProcA2C::run()/Timing: " << card2name_(card_) << endl;
1038 tm.Print();
1039 cout << " ---------------------------------------------------------- " << endl;
1040
1041 }
1042 catch (PException& exc) {
1043 cout << " BRProcA2C::run()/catched PException " << exc.Msg() << endl;
1044 setRC(3);
1045 return;
1046 }
1047 catch(...) {
1048 cout << " BRProcA2C::run()/catched unknown ... exception " << endl;
1049 setRC(4);
1050 return;
1051 }
1052 setRC(0);
1053 return;
1054}
1055
1056
1057//---------------------------------------------------------------------
1058// Classe thread de traitement 2 x 2 voies/frames (Apres BRProcA2C)
1059// !!!! OBSOLETE !!!!
1060//---------------------------------------------------------------------
1061/*!
1062 \class BRProcB4C
1063 \ingroup TAcq
1064
1065 \brief Deprecated processing class.
1066*/
1067
1068/* --Methode-- */
1069BRProcB4C::BRProcB4C(RAcqMemZoneMgr& mem1, RAcqMemZoneMgr& mem2, string& path,
1070 bool fgraw, uint_4 nmean, uint_4 nmax, bool fgnotrl)
1071 : memgr1(mem1), memgr2(mem2)
1072{
1073 fgraw_ = fgraw;
1074 nmax_ = nmax;
1075 nmean_ = nmean;
1076 if (fgraw_) cout << " BRProcB4C::BRProcB4C() - constructeur RAW data - NMean= " << nmean_ << endl;
1077 else cout << " BRProcB4C::BRProcB4C() - constructeur FFT data - NMean= " << nmean_ << endl;
1078 stop_ = false;
1079 path_ = path;
1080 fgnotrl_ = fgnotrl;
1081}
1082
1083/* --Methode-- */
1084void BRProcB4C::Stop()
1085{
1086 stop_=true;
1087 // cout <<" BRProcB4C::Stop ... > STOP " << endl;
1088}
1089
1090
1091/* --Methode-- */
1092void BRProcB4C::run()
1093{
1094 setRC(1);
1095 try {
1096 Timer tm("BRProcB4C", false);
1097 TimeStamp ts;
1098 BRPaqChecker pcheck1(!fgnotrl_); // Verification/comptage des paquets
1099 BRPaqChecker pcheck2(!fgnotrl_); // Verification/comptage des paquets
1100
1101 size_t totnbytesout = 0;
1102 size_t totnbytesproc = 0;
1103
1104 cout << " BRProcB4C::run() - Starting " << ts << " NMaxMemZones=" << nmax_
1105 << " NMean=" << nmean_ << endl;
1106 cout << " BRProcB4C::run()... - Output Data Path: " << path_ << endl;
1107
1108 uint_4 paqsz = memgr1.PaqSize();
1109 uint_4 procpaqsz = memgr1.ProcPaqSize();
1110 if ((paqsz != memgr2.PaqSize())||(procpaqsz!= memgr2.ProcPaqSize())) {
1111 cout << "BRProcB4C::run()/ERROR : different paquet size -> stop \n ...(PaqSz1="
1112 << paqsz << " Sz2=" << memgr2.PaqSize() << " ProcPaqSz1="
1113 << procpaqsz << " Sz2=" << memgr2.ProcPaqSize() << " )" << endl;
1114 setRC(9);
1115 return;
1116 }
1117
1118 TVector< complex<r_4> > cfour; // composant TF
1119 BRPaquet pq(NULL, NULL, paqsz);
1120/*
1121 TVector<r_4> vx(pq.DataSize()/2);
1122 vx = (r_4)(0.);
1123 FFTPackServer ffts;
1124 ffts.FFTForward(vx, cfour);
1125
1126 TVector< complex<r_4> > visiV13( cfour.Size() );
1127 TVector< complex<r_4> > visiV14( cfour.Size() );
1128 TVector< complex<r_4> > visiV23( cfour.Size() );
1129 TVector< complex<r_4> > visiV24( cfour.Size() );
1130*/
1131 int szfour = pq.DataSize()/2/2+1;
1132// int szfour = (paqsz-40)/2+1;
1133 TVector< complex<r_4> > visiV13( szfour );
1134 TVector< complex<r_4> > visiV14( szfour );
1135 TVector< complex<r_4> > visiV23( szfour );
1136 TVector< complex<r_4> > visiV24( szfour );
1137 // cout << " *DBG*AAAAA ---- Vectors OK" << endl;
1138 cout << " *DBG*BRProcB4C PaqSz=" << paqsz << " ProcPaqSize=" << procpaqsz
1139 << " procpaqsz/2=" << procpaqsz/2 << " cfour.Size()=" << szfour
1140 << " *8=" << szfour*8 << endl;
1141
1142 DataTable dt;
1143 dt.AddLongColumn("fc1");
1144 dt.AddLongColumn("tt1");
1145 dt.AddLongColumn("fc2");
1146 dt.AddLongColumn("tt2");
1147 DataTableRow dtr = dt.EmptyRow();
1148
1149 uint_4 nzm = 0;
1150 uint_4 totnoksfc = 0;
1151 uint_4 totnokpaq = 0;
1152 uint_4 totnpaq = 0;
1153 uint_4 ifile = 0;
1154
1155 uint_4 curfc=0;
1156 uint_8 curtt=0;
1157 uint_4 curfc2=0;
1158 uint_8 curtt2=0;
1159 uint_8 firsttt=0;
1160 uint_8 firsttt2=0;
1161 bool fgfirst=true;
1162 for (uint_4 kmz=0; kmz<nmax_; kmz++) {
1163 uint_4 noksfc = 0;
1164 uint_4 nokpaq = 0;
1165 if (stop_) break;
1166 // cout << " *DBG*BBBBB" << kmz << endl;
1167
1168 int mid1 = memgr1.FindMemZoneId(MemZA_ProcB);
1169 Byte* buff1 = memgr1.GetMemZone(mid1);
1170 if (buff1 == NULL) {
1171 cout << " BRProcB4C::run()/ERROR memgr.GetMemZone(" << mid1 << ") -> NULL" << endl;
1172 break;
1173 }
1174 Byte* procbuff1 = memgr1.GetProcMemZone(mid1);
1175 if (procbuff1 == NULL) {
1176 cout << " BRProcB4C::run()/ERROR memgr.GetProcMemZone(" << mid1 << ") -> NULL" << endl;
1177 break;
1178 }
1179 int mid2 = memgr2.FindMemZoneId(MemZA_ProcB);
1180 Byte* buff2 = memgr2.GetMemZone(mid2);
1181 if (buff1 == NULL) {
1182 cout << " BRProcB4C::run()/ERROR memgr.GetMemZone(" << mid2 << ") -> NULL" << endl;
1183 break;
1184 }
1185 Byte* procbuff2 = memgr2.GetProcMemZone(mid2);
1186 if (procbuff2 == NULL) {
1187 cout << " BRProcB4C::run()/ERROR memgr.GetProcMemZone(" << mid2 << ") -> NULL" << endl;
1188 break;
1189 }
1190 uint_4 i1,i2;
1191 i1=i2=0;
1192// cout << " *DBG*CCCCCC " << kmz << " memgr1.NbPaquets() =" << memgr1.NbPaquets() << endl;
1193 while((i1<memgr1.NbPaquets())&&(i2<memgr2.NbPaquets())) {
1194 BRPaquet paq1(NULL, buff1+i1*paqsz, paqsz);
1195 BRPaquet paq2(NULL, buff2+i2*paqsz, paqsz);
1196 totnpaq++;
1197// cout << " DBG["<<kmz<<"] i1,i2=" << i1 <<","<<i2<<" FC1,FC2=" <<paq1.FrameCounter()
1198//<<","<<paq2.FrameCounter()<<endl;
1199 // on ne traite que les paquets OK
1200 if (!pcheck1.Check(paq1)) { i1++; continue; }
1201 if (!pcheck2.Check(paq2)) { i2++; continue; }
1202 nokpaq++;
1203 if (paq1.FrameCounter()<paq2.FrameCounter()) { i1++; continue; }
1204 if (paq2.FrameCounter()<paq1.FrameCounter()) { i2++; continue; }
1205// cout << " DBG["<<kmz<<"]OKOK i1,i2=" << i1 <<","<<i2<<" FC1,FC2=" <<paq1.FrameCounter()
1206// <<","<<paq2.FrameCounter()<<endl;
1207
1208 if ((i1>=memgr1.NbPaquets())||(i2>=memgr1.NbPaquets())) {
1209 cout << " *BUG*BUG i1=" << i1 << " i2=" << i2 << endl;
1210 break;
1211 }
1212 // Les deux framecounters sont identiques ...
1213 noksfc++;
1214 curfc=paq1.FrameCounter();
1215 curfc2=paq2.FrameCounter();
1216 if (fgfirst) {
1217 firsttt=paq1.TimeTag(); firsttt2=paq2.TimeTag();
1218 cout << " BRProcB4C()/Info First FC="<<curfc<<" , "<<curfc2<<" -> TT="
1219 << firsttt<<" , "<<firsttt2 <<endl;
1220 fgfirst=false;
1221 }
1222 curtt=paq1.TimeTag()-firsttt;
1223 curtt2=paq2.TimeTag()-firsttt2;
1224 dtr[0]=curfc; dtr[1]=curtt;
1225 dtr[2]=curfc2; dtr[3]=curtt2;
1226 dt.AddRow(dtr);
1227
1228 complex<r_4>* zp1 = (complex<r_4>*)(procbuff1+i1*procpaqsz);
1229 complex<r_4>* zp2 = (complex<r_4>*)(procbuff1+i1*procpaqsz+procpaqsz/2);
1230 complex<r_4>* zp3 = (complex<r_4>*)(procbuff2+i2*procpaqsz);
1231 complex<r_4>* zp4 = (complex<r_4>*)(procbuff2+i2*procpaqsz+procpaqsz/2);
1232 for(sa_size_t j=0; j<visiV13.Size(); j++) {
1233 visiV13(j)+=zp1[j]*conj(zp3[j]);
1234 visiV14(j)+=zp1[j]*conj(zp4[j]);
1235 visiV23(j)+=zp2[j]*conj(zp3[j]);
1236 visiV24(j)+=zp2[j]*conj(zp4[j]);
1237 }
1238 if (nzm==0) {
1239 visiV13.Info()["StartFC"] = curfc;
1240 visiV14.Info()["StartFC"] = curfc;
1241 visiV23.Info()["StartFC"] = curfc;
1242 visiV24.Info()["StartFC"] = curfc;
1243 visiV13.Info()["StartTT"] = curtt;
1244 visiV14.Info()["StartTT"] = curtt;
1245 visiV23.Info()["StartTT"] = curtt;
1246 visiV24.Info()["StartTT"] = curtt;
1247 }
1248 nzm++; i1++; i2++;
1249 totnbytesproc += 2*paq1.DataSize();
1250 } // Fin de boucle sur les paquets d'une zone
1251 memgr1.FreeMemZone(mid1, MemZS_ProcB);
1252 memgr2.FreeMemZone(mid2, MemZS_ProcB);
1253
1254 if ((nzm >= nmean_) || ((kmz==(nmax_-1))&&(nzm>1))) {
1255 visiV13 /= complex<r_4>((r_4)nzm, 0.);
1256 visiV14 /= complex<r_4>((r_4)nzm, 0.);
1257 visiV23 /= complex<r_4>((r_4)nzm, 0.);
1258 visiV24 /= complex<r_4>((r_4)nzm, 0.);
1259 visiV13.Info()["NPaqMoy"] = nzm;
1260 visiV14.Info()["NPaqMoy"] = nzm;
1261 visiV23.Info()["NPaqMoy"] = nzm;
1262 visiV24.Info()["NPaqMoy"] = nzm;
1263 visiV13.Info()["EndFC"] = curfc;
1264 visiV14.Info()["EndFC"] = curfc;
1265 visiV23.Info()["EndFC"] = curfc;
1266 visiV24.Info()["EndFC"] = curfc;
1267 visiV13.Info()["EndTT"] = curtt;
1268 visiV14.Info()["EndTT"] = curtt;
1269 visiV23.Info()["EndTT"] = curtt;
1270 visiV24.Info()["EndTT"] = curtt;
1271 char fname[512];
1272 {
1273 sprintf(fname,"%s_%d.ppf",path_.c_str(),(int)ifile);
1274 POutPersist po(fname);
1275 po << PPFNameTag("visiV13") << visiV13;
1276 po << PPFNameTag("visiV14") << visiV14;
1277 po << PPFNameTag("visiV23") << visiV23;
1278 po << PPFNameTag("visiV24") << visiV24;
1279 }
1280 visiV13 = complex<r_4>(0., 0.);
1281 visiV14 = complex<r_4>(0., 0.);
1282 visiV23 = complex<r_4>(0., 0.);
1283 visiV24 = complex<r_4>(0., 0.);
1284 nzm = 0; ifile++;
1285// ts.SetNow();
1286// filog << ts << " : proc file " << fname << endl;
1287 cout << " BRProcB4C::run() created file " << fname << endl;
1288 }
1289 double okfrac = (nokpaq>1)?((double)noksfc/(double)nokpaq*100.):0.;
1290 cout << "BRProcB4C ["<<kmz<<"] NOKPaq=" << nokpaq << " NSameFC=" << noksfc
1291 << " (" << okfrac << " %)" << endl;
1292 totnokpaq += nokpaq;
1293 totnoksfc += noksfc;
1294 } // Fin de boucle sur les zones a traiter
1295 cout << " ------------------ BRProcB4C::run() END ----------------- " << endl;
1296 {
1297 dt.Info()["FirstTT1"]=firsttt;
1298 dt.Info()["FirstTT2"]=firsttt2;
1299 cout << dt;
1300 char fname[512];
1301 sprintf(fname,"%s_fctt.ppf",path_.c_str());
1302 POutPersist po(fname);
1303 po << PPFNameTag("ttfc") << dt;
1304 cout << " BRProcB4C::run() created TimeTag/FrameCounter file " << fname << endl;
1305 }
1306 ts.SetNow();
1307 tm.SplitQ();
1308 cout << " TotalProc= " << totnbytesproc/(1024*1024) << " MBytes, rate= "
1309 << (double)(totnbytesproc)/1024./tm.PartialElapsedTimems() << " MB/s" << endl;
1310 double totokfrac = (totnokpaq>1)?((double)totnoksfc/(double)totnokpaq*100.):0.;
1311 cout << " NOkPaq1,2=" << totnokpaq << " /TotNPaq=" << totnpaq << " TotNSameFC="
1312 << totnoksfc << " (" << totokfrac << " %)" << endl;
1313// cout << pcheck1;
1314// cout << pcheck2;
1315 cout << " BRProcB4C::run()/Timing: \n";
1316 tm.Print();
1317 cout << " ---------------------------------------------------------- " << endl;
1318}
1319 catch (PException& exc) {
1320 cout << " BRProcB4C::run()/catched PException " << exc.Msg() << endl;
1321 setRC(3);
1322 return;
1323 }
1324 catch(...) {
1325 cout << " BRProcB4C::run()/catched unknown ... exception " << endl;
1326 setRC(4);
1327 return;
1328 }
1329 setRC(0);
1330 return;
1331}
1332
1333
Note: See TracBrowser for help on using the repository browser.