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

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

1/ La description s'applique a la revision precedente qui a ete committee
avec un commentaire errone
2/ Mise en place des methodes et mecanismes permettant de traiter et transmettre correctement
le temps d'acquisition et le temps des paquets depuis la lecture fits aux fichiers de
traitement
3/ Ameliorations diverses ds BRMeanSpecCalculator et specmfib.cc

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