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

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

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