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

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

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