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

Last change on this file since 3703 was 3702, checked in by ansari, 16 years ago

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