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

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

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