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

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

Mise a jour et ajout de fichier pour taritement multifibres apres

prise de donnees de Nov2009 a Pittsburgh

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