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

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

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