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

Last change on this file since 3916 was 3905, checked in by ansari, 15 years ago

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