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

Last change on this file since 3693 was 3692, checked in by ansari, 16 years ago

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