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

Last change on this file since 3699 was 3698, checked in by ansari, 16 years ago

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