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

Last change on this file since 3956 was 3946, checked in by campagne, 15 years ago

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