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

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

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