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-- */
|
---|
33 | BRVisibilityCalculator::BRVisibilityCalculator(RAcqMemZoneMgr& memgr, string outpath, uint_4 nmean,
|
---|
34 | uint_4 freq1, uint_4 freq2, uint_4 nbfreq, size_t nthr)
|
---|
35 | : BRBaseProcessor(memgr), paralex_(*this, nthr), nparthr_(nthr), outpath_(outpath), nmean_(nmean),
|
---|
36 | numfreq1_(freq1), numfreq2_(freq2), nbinfreq_(nbfreq), 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 | fgallfibok=NULL;
|
---|
79 | fgcktt_=false;
|
---|
80 | }
|
---|
81 |
|
---|
82 | /* --Methode-- */
|
---|
83 | BRVisibilityCalculator::~BRVisibilityCalculator()
|
---|
84 | {
|
---|
85 | cout << " BRVisibilityCalculator - Visibility Datatable : " << endl;
|
---|
86 | cout << visdt_;
|
---|
87 | POutPersist po(outpath_+"visdt.ppf");
|
---|
88 | po << visdt_;
|
---|
89 | POutPersist poc(outpath_+"chanum.ppf");
|
---|
90 | poc << chanum_;
|
---|
91 |
|
---|
92 | if (fgcktt_) {
|
---|
93 | cout << " BRVisibilityCalculator - Check TimeTag done: TotNPaqProc= " << totnbpaq_ << endl;
|
---|
94 | for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
|
---|
95 | cout << " BRTTCheck-Fiber[" << fib << "] NBadTT=" << vbadtt_[fib] << " NDiffTT>5="
|
---|
96 | << vndiff5tt_[fib] << " NotSameTT=" << vnsamett_[fib] << endl;
|
---|
97 | }
|
---|
98 | POutPersist pott(outpath_+"ttfcmtx.ppf");
|
---|
99 | pott << PPFNameTag("FC") << fcmtx_;
|
---|
100 | pott << PPFNameTag("TT") << ttmtx_;
|
---|
101 | }
|
---|
102 | delete[] xnt_;
|
---|
103 | }
|
---|
104 |
|
---|
105 | /* --Methode-- */
|
---|
106 | void BRVisibilityCalculator::run()
|
---|
107 | {
|
---|
108 | if (nparthr_ < 2) return BRBaseProcessor::run();
|
---|
109 | // Execution multithread parallele
|
---|
110 | setRC(1);
|
---|
111 | int rc=0;
|
---|
112 | try {
|
---|
113 | cout << " BRVisibilityCalculator::run() - Starting " << " NFibers=" << memgr_.NbFibres()
|
---|
114 | << " NChan=" << 2*memgr_.NbFibres() << endl;
|
---|
115 | if ((nmean_%memgr_.NbPaquets())!=0) {
|
---|
116 | uint_4 mnmean = (nmean_/memgr_.NbPaquets()+1)*memgr_.NbPaquets();
|
---|
117 | cout << " BRVisibilityCalculator::run()/Info changing nmean=" << nmean_ << " to multiple of"
|
---|
118 | << " memgr_.NbPaquets() -> " << mnmean << endl;
|
---|
119 | nmean_=mnmean;
|
---|
120 | }
|
---|
121 | paralex_.SetParallelTask(*this);
|
---|
122 | cout << " BRVisibilityCalculator::run()/Info : starting ParallelExecutor with nThreads="
|
---|
123 | << paralex_.nThreads() << " ... " << endl;
|
---|
124 | paralex_.start();
|
---|
125 |
|
---|
126 | fgallfibok = new bool[memgr_.NbPaquets()];
|
---|
127 |
|
---|
128 | size_t paqsz=memgr_.PaqSize();
|
---|
129 | bool fgrun=true;
|
---|
130 | while (fgrun) {
|
---|
131 | if (stop_) break;
|
---|
132 | if (memgr_.GetRunState() == MemZR_Stopped) break;
|
---|
133 | int mid = memgr_.FindMemZoneId(mmact_); // (MemZA_ProcA);
|
---|
134 | Byte* buffg = memgr_.GetMemZone(mid);
|
---|
135 | if (buffg == NULL) {
|
---|
136 | cout << "BRVisibilityCalculator::run()/ERROR memgr.GetMemZone(" << mid << ") -> NULL" << endl;
|
---|
137 | setRC(7); fgrun=false;
|
---|
138 | break;
|
---|
139 | }
|
---|
140 | for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
|
---|
141 | fbuff_[fib] = memgr_.GetMemZone(mid,fib);
|
---|
142 | if (fbuff_[fib] == NULL) { // cela ne devrait pas arriver
|
---|
143 | cout << "BRVisibilityCalculator::run()/ERROR memgr.GetMemZone(" << mid << "," << fib << ") -> NULL" << endl;
|
---|
144 | setRC(9); fgrun=false;
|
---|
145 | break;
|
---|
146 | }
|
---|
147 | }
|
---|
148 |
|
---|
149 | if (totnbpaq_%nmean_ == 0) {
|
---|
150 | if (totnbpaq_ > 0) {
|
---|
151 | moyfc_/=nmean_;
|
---|
152 | moytt_/=nmean_;
|
---|
153 | vismtx_.Info()["MeanFC"] = moyfc_;
|
---|
154 | vismtx_.Info()["MeanTT"] = moytt_;
|
---|
155 | vismtx_.Info()["NPAQSUM"] = nmean_;
|
---|
156 |
|
---|
157 | // ATTENTION : Matrice visibilites non moyennee
|
---|
158 | char nfile[32];
|
---|
159 | sprintf(nfile,"vismtx%d.ppf",numfile_);
|
---|
160 | string flnm=outpath_+nfile;
|
---|
161 | POutPersist po(flnm);
|
---|
162 | po << vismtx_;
|
---|
163 | cout << numfile_ << "-BRVisibilityCalculator::run() NPaqProc="
|
---|
164 | << totnbpaq_ << " -> Visibility Matrix in " << flnm << endl;
|
---|
165 | FillVisibTable(moyfc_, moytt_);
|
---|
166 | numfile_++;
|
---|
167 | }
|
---|
168 | vismtx_ = complex<r_4>((r_4)0.,(r_4)0.);
|
---|
169 | moyfc_=moytt_=0.;
|
---|
170 | }
|
---|
171 |
|
---|
172 | for(size_t jp=0; jp<memgr_.NbPaquets(); jp++) { // boucle sur les paquets d'une zone
|
---|
173 | fgallfibok[jp]=fgokallfibers_=true;
|
---|
174 | for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
|
---|
175 | vpaq_[fib].Set(fbuff_[fib]+jp*paqsz);
|
---|
176 | vfgok_[fib] = vpchk_[fib].Check(vpaq_[fib],curfc_[fib]);
|
---|
177 | if (!vfgok_[fib]) fgallfibok[jp]=fgokallfibers_=false;
|
---|
178 | }
|
---|
179 | if (fgokallfibers_) {
|
---|
180 | if (totprocnpaq_==0) {
|
---|
181 | for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
|
---|
182 | fcfirst_[fib]=curfc_[fib];
|
---|
183 | ttfirst_[fib]=vpaq_[fib].TimeTag();
|
---|
184 | }
|
---|
185 | }
|
---|
186 | totprocnpaq_++;
|
---|
187 | moyfc_ += curfc_[0];
|
---|
188 | moytt_ += (vpaq_[0].TimeTag()-ttfirst_[0]);
|
---|
189 | if (fgcktt_) CheckTimeTag();
|
---|
190 | totnbpaq_++;
|
---|
191 | }
|
---|
192 | } // Fin de boucle sur les paquets
|
---|
193 |
|
---|
194 | // Execution parallele pour calcul des visibilites par bandes de frequence
|
---|
195 | int rcpex=paralex_.execute();
|
---|
196 | if (rcpex!=0) cout << " BRVisibilityCalculator::run() / Error Rc[paralex_.execute()]=" << rcpex << endl;
|
---|
197 |
|
---|
198 | memgr_.FreeMemZone(mid, mmsta_); // (MemZS_ProcA);
|
---|
199 | } // Fin de boucle sur les zones a traiter
|
---|
200 | //------------------------------------
|
---|
201 | cout << " --------- END BRVisibilityCalculator::run() , TotNbProcPaq=" << totprocnpaq_ << endl;
|
---|
202 | /*
|
---|
203 | for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) vpchk_[fib].Print();
|
---|
204 | cout << " ------------------------------------ " << endl;
|
---|
205 | */
|
---|
206 | delete[] fgallfibok;
|
---|
207 | }
|
---|
208 | catch (std::exception& exc) {
|
---|
209 | cout << " BRVisibilityCalculator::run()/catched std::exception " << exc.what() << endl;
|
---|
210 | setRC(98);
|
---|
211 | return;
|
---|
212 | }
|
---|
213 | catch(...) {
|
---|
214 | cout << " BRVisibilityCalculator::run()/catched unknown ... exception " << endl;
|
---|
215 | setRC(99);
|
---|
216 | return;
|
---|
217 | }
|
---|
218 |
|
---|
219 | }
|
---|
220 |
|
---|
221 | /* --Methode-- */
|
---|
222 | int BRVisibilityCalculator::ActivateTimeTagCheck(uint_8 maxnpaq)
|
---|
223 | {
|
---|
224 | mindeltatt_=memgr_.PaqSize()/2;
|
---|
225 | if (mindeltatt_<1) mindeltatt_=1;
|
---|
226 | fcmtx_.SetSize(memgr_.NbFibres(), maxnpaq);
|
---|
227 | ttmtx_.SetSize(memgr_.NbFibres(), maxnpaq);
|
---|
228 | vlasttt_.resize(memgr_.NbFibres(), 0);
|
---|
229 | vbadtt_.resize(memgr_.NbFibres(), 0);
|
---|
230 | vnsamett_.resize(memgr_.NbFibres(), 0);
|
---|
231 | vndiff5tt_.resize(memgr_.NbFibres(), 0);
|
---|
232 |
|
---|
233 | fgcktt_=true;
|
---|
234 | cout << " BRVisibilityCalculator::ActivateTimeTagCheck() - TT/Fc matrix NCols=" << maxnpaq
|
---|
235 | << " MinDeltaTT=" << mindeltatt_ << endl;
|
---|
236 |
|
---|
237 | return 0;
|
---|
238 | }
|
---|
239 |
|
---|
240 | /* --Methode-- */
|
---|
241 | int BRVisibilityCalculator::Process()
|
---|
242 | {
|
---|
243 |
|
---|
244 | for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
|
---|
245 | vpdata_[2*fib] = vpaq_[fib].Data1C();
|
---|
246 | vpdata_[2*fib+1] = vpaq_[fib].Data2C();
|
---|
247 | }
|
---|
248 |
|
---|
249 | if (totnbpaq_%nmean_ == 0) {
|
---|
250 | if (totnbpaq_ > 0) {
|
---|
251 | moyfc_/=nmean_;
|
---|
252 | moytt_/=nmean_;
|
---|
253 | vismtx_.Info()["MeanFC"] = moyfc_;
|
---|
254 | vismtx_.Info()["MeanTT"] = moytt_;
|
---|
255 | vismtx_.Info()["NPAQSUM"] = nmean_;
|
---|
256 |
|
---|
257 | // ATTENTION : Matrice visibilites non moyennee
|
---|
258 | char nfile[32];
|
---|
259 | sprintf(nfile,"vismtx%d.ppf",numfile_);
|
---|
260 | string flnm=outpath_+nfile;
|
---|
261 | POutPersist po(flnm);
|
---|
262 | po << vismtx_;
|
---|
263 | cout << numfile_ << "-BRVisibilityCalculator::Process() NPaqProc="
|
---|
264 | << totnbpaq_ << " -> Visibility Matrix in " << flnm << endl;
|
---|
265 | FillVisibTable(moyfc_, moytt_);
|
---|
266 | numfile_++;
|
---|
267 | }
|
---|
268 | vismtx_ = complex<r_4>((r_4)0.,(r_4)0.);
|
---|
269 | moyfc_=moytt_=0.;
|
---|
270 | }
|
---|
271 |
|
---|
272 | sa_size_t k=0;
|
---|
273 | for(size_t i=0; i<vpdata_.size(); i++) {
|
---|
274 | for(size_t j=i; j<vpdata_.size(); j++) {
|
---|
275 | TVector< complex<r_4> > vis = vismtx_.Row(k); k++;
|
---|
276 | for(sa_size_t f=1; f<vis.Size(); f++) {
|
---|
277 | vis(f) += complex<r_4>((r_4)vpdata_[i][f].realB(), (r_4)vpdata_[i][f].imagB()) *
|
---|
278 | complex<r_4>((r_4)vpdata_[j][f].realB(), -(r_4)vpdata_[j][f].imagB());
|
---|
279 | }
|
---|
280 | }
|
---|
281 | }
|
---|
282 |
|
---|
283 | moyfc_ += curfc_[0];
|
---|
284 | moytt_ += (vpaq_[0].TimeTag()-ttfirst_[0]);
|
---|
285 | if (fgcktt_) CheckTimeTag();
|
---|
286 | totnbpaq_++;
|
---|
287 | return 0;
|
---|
288 | }
|
---|
289 |
|
---|
290 | /* --Methode-- */
|
---|
291 | int BRVisibilityCalculator::execute(int tid)
|
---|
292 | {
|
---|
293 | vector<TwoByteComplex*> pvpdata(2*memgr_.NbFibres());
|
---|
294 | size_t paqsz=memgr_.PaqSize();
|
---|
295 | BRPaquet ppaq(paqsz);
|
---|
296 |
|
---|
297 | sa_size_t fdelt = vismtx_.NCols()/nparthr_;
|
---|
298 | sa_size_t fdeb = tid*fdelt;
|
---|
299 | sa_size_t ffin = (tid+1)*fdelt;
|
---|
300 |
|
---|
301 | if (fdeb<1) fdeb=1;
|
---|
302 | if ((ffin>vismtx_.NCols())||(tid==(nparthr_-1))) ffin=vismtx_.NCols();
|
---|
303 |
|
---|
304 | for(size_t jp=0; jp<memgr_.NbPaquets(); jp++) { // boucle sur les paquets d'une zone
|
---|
305 | if (!fgallfibok[jp]) continue;
|
---|
306 | for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
|
---|
307 | ppaq.Set(fbuff_[fib]+jp*paqsz);
|
---|
308 | pvpdata[2*fib] = ppaq.Data1C();
|
---|
309 | pvpdata[2*fib+1] = ppaq.Data2C();
|
---|
310 | }
|
---|
311 | sa_size_t k=0;
|
---|
312 | for(size_t i=0; i<vpdata_.size(); i++) {
|
---|
313 | for(size_t j=i; j<vpdata_.size(); j++) {
|
---|
314 | TVector< complex<r_4> > vis = vismtx_.Row(k); k++;
|
---|
315 | for(sa_size_t f=fdeb; f<ffin; f++) {
|
---|
316 | vis(f) += complex<r_4>((r_4)pvpdata[i][f].realB(), (r_4)pvpdata[i][f].imagB()) *
|
---|
317 | complex<r_4>((r_4)pvpdata[j][f].realB(), -(r_4)pvpdata[j][f].imagB());
|
---|
318 | }
|
---|
319 | }
|
---|
320 | }
|
---|
321 |
|
---|
322 | } // Fin de boucle sur les paquets
|
---|
323 |
|
---|
324 | return 0;
|
---|
325 | }
|
---|
326 |
|
---|
327 | /* --Methode-- */
|
---|
328 | int BRVisibilityCalculator::FillVisibTable(double fcm, double ttm)
|
---|
329 | {
|
---|
330 | xnt_[0]=fcm; xnt_[1]=ttm/1.25e8;
|
---|
331 |
|
---|
332 | if (djf_<2) {
|
---|
333 | for(sa_size_t rv=0; rv<vismtx_.NRows(); rv++) {
|
---|
334 | for(sa_size_t jf=jf1_; jf<jf2_; jf++) {
|
---|
335 | xnt_[2]=jf;
|
---|
336 | xnt_[3]=chanum_(rv);
|
---|
337 | xnt_[4]=vismtx_(rv,jf).real()/(r_4)(nmean_);
|
---|
338 | xnt_[5]=vismtx_(rv,jf).imag()/(r_4)(nmean_);
|
---|
339 | visdt_.AddRow(xnt_);
|
---|
340 | }
|
---|
341 | }
|
---|
342 | }
|
---|
343 | else {
|
---|
344 | for(sa_size_t rv=0; rv<vismtx_.NRows(); rv++) {
|
---|
345 | for(sa_size_t jf=jf1_; jf<jf2_; jf+=djf_) {
|
---|
346 | r_4 moyreal=0.;
|
---|
347 | r_4 moyimag=0.;
|
---|
348 | sa_size_t jjfmx=jf+djf_;
|
---|
349 | if (jjfmx > vismtx_.NCols()) jjfmx=vismtx_.NCols();
|
---|
350 | for(sa_size_t jjf=jf; jjf<jjfmx; jjf++) {
|
---|
351 | moyreal+=vismtx_(rv,jjf).real();
|
---|
352 | moyimag+=vismtx_(rv,jjf).imag();
|
---|
353 | }
|
---|
354 | xnt_[2]=jf+djf_/2;
|
---|
355 | xnt_[3]=chanum_(rv);
|
---|
356 | xnt_[4]=moyreal/(r_4)(nmean_*djf_);
|
---|
357 | xnt_[5]=moyimag/(r_4)(nmean_*djf_);
|
---|
358 | visdt_.AddRow(xnt_);
|
---|
359 | }
|
---|
360 | }
|
---|
361 | }
|
---|
362 | return 0;
|
---|
363 | }
|
---|
364 |
|
---|
365 | /* --Methode-- */
|
---|
366 | int BRVisibilityCalculator::CheckTimeTag()
|
---|
367 | {
|
---|
368 | if (totnbpaq_==0) {
|
---|
369 | for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
|
---|
370 | vlasttt_[fib]=ttfirst_[fib];
|
---|
371 | if (ttmtx_.NCols()>0) {
|
---|
372 | fcmtx_(fib,totnbpaq_) = curfc_[fib];
|
---|
373 | ttmtx_(fib,totnbpaq_) = vlasttt_[fib];
|
---|
374 | }
|
---|
375 | }
|
---|
376 | return 0;
|
---|
377 | }
|
---|
378 | for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
|
---|
379 | int_8 ld = (int_8)vpaq_[fib].TimeTag()-(int_8)vlasttt_[fib];
|
---|
380 | int_8 fd = (int_8)vpaq_[fib].TimeTag()-(int_8)ttfirst_[fib]-(int_8)vpaq_[0].TimeTag()+(int_8)ttfirst_[0];
|
---|
381 | /* if ( (ld < mindeltatt_) || (fd<-5) || (fd>5)) { vbadtt_[fib]++; vnsamett_[fib]++; }
|
---|
382 | else {
|
---|
383 | if (fd!=0) vnsamett_[fib]++;
|
---|
384 | }
|
---|
385 | */
|
---|
386 | if (ld < mindeltatt_) vbadtt_[fib]++;
|
---|
387 | else {
|
---|
388 | if (fd != 0) vnsamett_[fib]++;
|
---|
389 | if ((fd<-5)||(fd>5)) vndiff5tt_[fib]++;
|
---|
390 | }
|
---|
391 | vlasttt_[fib]=vpaq_[fib].TimeTag();
|
---|
392 | if (totnbpaq_<ttmtx_.NCols()) {
|
---|
393 | fcmtx_(fib,totnbpaq_) = curfc_[fib];
|
---|
394 | ttmtx_(fib,totnbpaq_) = vlasttt_[fib];
|
---|
395 | }
|
---|
396 | }
|
---|
397 | return 0;
|
---|
398 | }
|
---|
399 |
|
---|
400 |
|
---|
401 | //---------------------------------------------------------------------
|
---|
402 | // Classe de traitement simple - calcul de spectres moyennes / voie
|
---|
403 | //---------------------------------------------------------------------
|
---|
404 | /* --Methode-- */
|
---|
405 | BRMeanSpecCalculator::BRMeanSpecCalculator(RAcqMemZoneMgr& memgr, string outpath, uint_4 nmean)
|
---|
406 | : BRBaseProcessor(memgr), outpath_(outpath), nmean_(nmean)
|
---|
407 | {
|
---|
408 | BRPaquet paq(memgr_.PaqSize());
|
---|
409 | mspecmtx_.SetSize(2*memgr_.NbFibres(), paq.DataSize()/4);
|
---|
410 | numfile_=0;
|
---|
411 | totnbpaq_=0;
|
---|
412 | }
|
---|
413 |
|
---|
414 | /* --Methode-- */
|
---|
415 | BRMeanSpecCalculator::~BRMeanSpecCalculator()
|
---|
416 | {
|
---|
417 | }
|
---|
418 |
|
---|
419 |
|
---|
420 | /* --Methode-- */
|
---|
421 | int BRMeanSpecCalculator::Process()
|
---|
422 | {
|
---|
423 |
|
---|
424 | if (totnbpaq_%nmean_ == 0) {
|
---|
425 | if (totnbpaq_ > 0) {
|
---|
426 | mspecmtx_.Info()["NPAQSUM"] = nmean_;
|
---|
427 | mspecmtx_ /= (double)nmean_;
|
---|
428 | char nfile[32];
|
---|
429 | sprintf(nfile,"mspecmtx%d.ppf",numfile_);
|
---|
430 | string flnm=outpath_+nfile;
|
---|
431 | POutPersist po(flnm);
|
---|
432 | po << mspecmtx_;
|
---|
433 | cout << numfile_ << "-BRMeanSpecCalculator::Process() NPaqProc="
|
---|
434 | << totnbpaq_ << " -> Mean spectra Matrix in " << flnm << endl;
|
---|
435 | numfile_++;
|
---|
436 | }
|
---|
437 | mspecmtx_ = (r_8)(0.);
|
---|
438 | }
|
---|
439 |
|
---|
440 | sa_size_t k=0;
|
---|
441 | for(size_t i=0; i<(size_t)2*memgr_.NbFibres(); i++) {
|
---|
442 | TwoByteComplex* zp=vpaq_[i/2].Data1C();
|
---|
443 | if (i%2==1) zp=vpaq_[i/2].Data2C();
|
---|
444 | TVector< r_4 > spec = mspecmtx_.Row(k); k++;
|
---|
445 | for(sa_size_t f=1; f<spec.Size(); f++) {
|
---|
446 | spec(f) += zp[f].module2F();
|
---|
447 | }
|
---|
448 | }
|
---|
449 |
|
---|
450 | totnbpaq_++;
|
---|
451 | return 0;
|
---|
452 | }
|
---|
453 |
|
---|
454 |
|
---|
455 | //---------------------------------------------------------------
|
---|
456 | // Classe thread de traitement donnees ADC avec 2 voies par frame
|
---|
457 | //---------------------------------------------------------------
|
---|
458 |
|
---|
459 | // Mutex pour eviter le plantage du a FFTW qui ne semble pas thread-safe
|
---|
460 | static ZMutex* pmutfftw=NULL;
|
---|
461 |
|
---|
462 | /* --Methode-- */
|
---|
463 | BRProcA2C::BRProcA2C(RAcqMemZoneMgr& mem, string& path, bool fgraw, uint_4 nmean,
|
---|
464 | uint_4 nmax, bool fghist, uint_4 nfsmap, bool fgnotrl, int card)
|
---|
465 | : memgr(mem)
|
---|
466 | {
|
---|
467 | fgraw_ = fgraw;
|
---|
468 | nmax_ = nmax;
|
---|
469 | nmean_ = nmean;
|
---|
470 | if (fgraw_) cout << " BRProcA2C::BRProcA2C() - constructeur RAW data - NMean=" << nmean_ << endl;
|
---|
471 | else cout << " BRProcA2C::BRProcA2C() - constructeur FFT data - NMean=" << nmean_ << endl;
|
---|
472 | nfsmap_ = nfsmap;
|
---|
473 | stop_ = false;
|
---|
474 | path_ = path;
|
---|
475 | fgnotrl_ = fgnotrl;
|
---|
476 | fghist_ = fghist;
|
---|
477 | card_ = card;
|
---|
478 | if (pmutfftw==NULL) pmutfftw=new ZMutex;
|
---|
479 | }
|
---|
480 |
|
---|
481 | /* --Methode-- */
|
---|
482 | void BRProcA2C::Stop()
|
---|
483 | {
|
---|
484 | stop_=true;
|
---|
485 | // cout <<" BRProcA2C::Stop ... > STOP " << endl;
|
---|
486 | }
|
---|
487 |
|
---|
488 |
|
---|
489 | static inline r_4 Zmod2(complex<r_4> z)
|
---|
490 | { return (z.real()*z.real()+z.imag()*z.imag()); }
|
---|
491 |
|
---|
492 | static inline string card2name_(int card)
|
---|
493 | {
|
---|
494 | if (card==2) return " (Chan3,4) ";
|
---|
495 | else return " (Chan1,2) ";
|
---|
496 | }
|
---|
497 | /* --Methode-- */
|
---|
498 | void BRProcA2C::run()
|
---|
499 | {
|
---|
500 | setRC(1);
|
---|
501 | try {
|
---|
502 | Timer tm("BRProcA2C", false);
|
---|
503 | TimeStamp ts;
|
---|
504 | BRPaqChecker pcheck(!fgnotrl_); // Verification/comptage des paquets
|
---|
505 |
|
---|
506 | size_t totnbytesout = 0;
|
---|
507 | size_t totnbytesproc = 0;
|
---|
508 |
|
---|
509 | cout << " BRProcA2C::run() - Starting " << ts << " NMaxMemZones=" << nmax_
|
---|
510 | << " NMean=" << nmean_ << card2name_(card_) << endl;
|
---|
511 | cout << " BRProcA2C::run()... - Output Data Path: " << path_ << endl;
|
---|
512 | char fname[512];
|
---|
513 | // sprintf(fname,"%s/proc.log",path_.c_str());
|
---|
514 | // ofstream filog(fname);
|
---|
515 | // filog << " BRProcA2C::run() - starting log file " << ts << endl;
|
---|
516 | // filog << " ... NMaxMemZones=" << nmax_ << " NMean=" << nmean_ << " Step=" << step_ << endl;
|
---|
517 |
|
---|
518 | /*----DELETE NTuple
|
---|
519 | const char* nnames[8] = {"fcs","tts","s1","s2","s12","s12re","s12im","s12phi"};
|
---|
520 | NTuple nt(8, nnames);
|
---|
521 | double xnt[10];
|
---|
522 | uint_4 nmnt = 0;
|
---|
523 | double ms1,ms2,ms12,ms12re,ms12im,ms12phi;
|
---|
524 | ----*/
|
---|
525 | // Time sample (raw data) /FFT coeff histograms
|
---|
526 | Histo* ph1=NULL;
|
---|
527 | Histo* ph2=NULL;
|
---|
528 | if (fghist_) {
|
---|
529 | if (fgraw_) {
|
---|
530 | ph1 = new Histo(-0.5, 255.5, 256);
|
---|
531 | ph2 = new Histo(-0.5, 255.5, 256);
|
---|
532 | }
|
---|
533 | else {
|
---|
534 | ph1 = new Histo(-128.5, 128.5, 257);
|
---|
535 | ph2 = new Histo(-128.5, 128.5, 257);
|
---|
536 | }
|
---|
537 | }
|
---|
538 |
|
---|
539 | // Initialisation pour calcul FFT
|
---|
540 | TVector< complex<r_4> > cfour1; // composant TF
|
---|
541 | uint_4 paqsz = memgr.PaqSize();
|
---|
542 | uint_4 procpaqsz = memgr.ProcPaqSize();
|
---|
543 |
|
---|
544 |
|
---|
545 | BRPaquet pq(NULL, NULL, paqsz);
|
---|
546 | TVector<r_4> vx(pq.DataSize()/2);
|
---|
547 | int szfour = pq.DataSize()/2/2+1;
|
---|
548 | cfour1.SetSize(szfour);
|
---|
549 | /*
|
---|
550 | vx = (r_4)(0.);
|
---|
551 | FFTPackServer ffts;
|
---|
552 | ffts.FFTForward(vx, cfour1);
|
---|
553 | szfour = cfour1.Size();
|
---|
554 | */
|
---|
555 |
|
---|
556 | bool fgtimfreq = false; // true->cartes temps<>frequences
|
---|
557 | if (nfsmap_>0) fgtimfreq=true;
|
---|
558 |
|
---|
559 | TVector< complex<r_4> > cfour2(cfour1.Size());
|
---|
560 |
|
---|
561 | TVector<r_4> spectreV1(cfour1.Size());
|
---|
562 | TVector<r_4> spectreV2(cfour1.Size());
|
---|
563 | TVector<r_4> moyspecV1(cfour1.Size()); // Moyenne des Spectres
|
---|
564 | TVector<r_4> moyspecV2(cfour1.Size());
|
---|
565 | TVector<r_4> sigspecV1(cfour1.Size()); // Sigma des Spectres
|
---|
566 | TVector<r_4> sigspecV2(cfour1.Size());
|
---|
567 | TVector< complex<r_4> > visiV12( cfour1.Size() );
|
---|
568 |
|
---|
569 | TMatrix<r_4> timfreqV1, timfreqV2; // Cartes temps<>frequences
|
---|
570 | if (fgtimfreq) {
|
---|
571 | timfreqV1.SetSize(nmean_, spectreV1.Size()/nfsmap_);
|
---|
572 | timfreqV2.SetSize(nmean_, spectreV2.Size()/nfsmap_);
|
---|
573 | }
|
---|
574 | cout << " *DBG*BRProcA2C PaqSz=" << paqsz << " ProcPaqSize=" << procpaqsz
|
---|
575 | << " procpaqsz/2=" << procpaqsz/2 << " cfour1.Size()=" << cfour1.Size()
|
---|
576 | << " *8=" << cfour1.Size()*8 << endl;
|
---|
577 |
|
---|
578 | pmutfftw->lock();
|
---|
579 | fftwf_plan plan1 = fftwf_plan_dft_r2c_1d(vx.Size(), vx.Data(),
|
---|
580 | (fftwf_complex*)cfour1.Data(), FFTW_ESTIMATE);
|
---|
581 | fftwf_plan plan2 = fftwf_plan_dft_r2c_1d(vx.Size(), vx.Data(),
|
---|
582 | (fftwf_complex*)cfour2.Data(), FFTW_ESTIMATE);
|
---|
583 | pmutfftw->unlock();
|
---|
584 |
|
---|
585 | uint_4 ifile = 0;
|
---|
586 | uint_4 nzm = 0; // Nb de paquets moyennes pour le calcul de chaque spectre
|
---|
587 | uint_4 nmoyspec = 0; // Nb de spectres moyennes
|
---|
588 |
|
---|
589 | uint_4 curfc=0;
|
---|
590 | uint_8 curtt=0;
|
---|
591 | uint_8 firsttt=0;
|
---|
592 | bool fgfirst=true;
|
---|
593 | double moysig[2]={0.,0.};
|
---|
594 | double sigsig[2]={0.,0.};
|
---|
595 | uint_8 nbsig[2]={0,0};
|
---|
596 |
|
---|
597 | for (uint_4 kmz=0; kmz<nmax_; kmz++) {
|
---|
598 | if (stop_) break;
|
---|
599 | int mid = memgr.FindMemZoneId(MemZA_ProcA);
|
---|
600 | Byte* buff = memgr.GetMemZone(mid);
|
---|
601 | if (buff == NULL) {
|
---|
602 | cout << " BRProcA2C::run()/ERROR memgr.GetMemZone(" << mid << ") -> NULL" << endl;
|
---|
603 | break;
|
---|
604 | }
|
---|
605 | Byte* procbuff = memgr.GetProcMemZone(mid);
|
---|
606 | if (procbuff == NULL) {
|
---|
607 | cout << " BRProcA2C::run()/ERROR memgr.GetProcMemZone(" << mid << ") -> NULL" << endl;
|
---|
608 | break;
|
---|
609 | }
|
---|
610 | //---- DELETE nmnt=0; ms1=ms2=ms12=ms12re=ms12im=ms12phi=0.;
|
---|
611 | for(uint_4 i=0; i<memgr.NbPaquets(); i++) {
|
---|
612 | BRPaquet paq(NULL, buff+i*paqsz, paqsz);
|
---|
613 | if (!pcheck.Check(paq)) continue; // on ne traite que les paquets OK
|
---|
614 | if (fgfirst) { firsttt=paq.TimeTag(); fgfirst=false; }
|
---|
615 | curfc=paq.FrameCounter();
|
---|
616 | curtt=paq.TimeTag()-firsttt;
|
---|
617 | // Traitement voie 1
|
---|
618 | if (fghist_) {
|
---|
619 | for(sa_size_t j=0; j<vx.Size(); j++) {
|
---|
620 | r_4 vts=(fgraw_)?((r_4)(*(paq.Data1()+j))):((r_4)(*(paq.Data1S()+j)));
|
---|
621 | ph1->Add((r_8)vts);
|
---|
622 | moysig[0] += (double)vts;
|
---|
623 | sigsig[0] += ((double)vts)*((double)vts);
|
---|
624 | nbsig[0]++;
|
---|
625 | }
|
---|
626 | for(sa_size_t j=0; j<vx.Size(); j++) {
|
---|
627 | r_4 vts=(fgraw_)?((r_4)(*(paq.Data2()+j))):((r_4)(*(paq.Data2S()+j)));
|
---|
628 | ph2->Add((r_8)vts);
|
---|
629 | moysig[1] += (double)vts;
|
---|
630 | sigsig[1] += ((double)vts)*((double)vts);
|
---|
631 | nbsig[1]++;
|
---|
632 | }
|
---|
633 | }
|
---|
634 | if (fgraw_) {
|
---|
635 | for(sa_size_t j=0; j<vx.Size(); j++)
|
---|
636 | vx(j) = (r_4)(*(paq.Data1()+j))-127.5;
|
---|
637 | // fftwf_complex* coeff1 = (fftwf_complex*)(procbuff+i*procpaqsz);
|
---|
638 | fftwf_execute(plan1);
|
---|
639 | // Traitement voie 2
|
---|
640 | for(sa_size_t j=0; j<vx.Size(); j++)
|
---|
641 | vx(j) = (r_4)(*(paq.Data2()+j))-127.5;
|
---|
642 | fftwf_execute(plan2);
|
---|
643 | }
|
---|
644 | else {
|
---|
645 | for(sa_size_t j=1; j<cfour1.Size()-1; j++) {
|
---|
646 | cfour1(j) = complex<r_4>((r_4)paq.Data1C()[j].realB(), (r_4)paq.Data1C()[j].imagB());
|
---|
647 | cfour2(j) = complex<r_4>((r_4)paq.Data2C()[j].realB(), (r_4)paq.Data2C()[j].imagB());
|
---|
648 | }
|
---|
649 | cfour1(0) = complex<r_4>((r_4)paq.Data1C()[0].realB(), (r_4)0.);
|
---|
650 | cfour1(cfour1.Size()-1) = complex<r_4>((r_4)paq.Data1C()[0].imagB(), (r_4)0.);
|
---|
651 | cfour2(0) = complex<r_4>((r_4)paq.Data2C()[0].realB(), (r_4)0.);
|
---|
652 | cfour2(cfour2.Size()-1) = complex<r_4>((r_4)paq.Data2C()[0].imagB(), (r_4)0.);
|
---|
653 | }
|
---|
654 | for(sa_size_t j=0; j<spectreV1.Size(); j++)
|
---|
655 | spectreV1(j) += Zmod2(cfour1(j));
|
---|
656 | memcpy(procbuff+i*procpaqsz, cfour1.Data(), sizeof(complex<r_4>)*cfour1.Size());
|
---|
657 | if (fgtimfreq) { // Remplissage tableau temps-frequence
|
---|
658 | for(sa_size_t c=1; c<timfreqV1.NCols(); c++) {
|
---|
659 | for(sa_size_t j=c*nfsmap_; j<(c+1)*nfsmap_; j++)
|
---|
660 | timfreqV1(nzm, c) += Zmod2(cfour1(j));
|
---|
661 | }
|
---|
662 | }
|
---|
663 | for(sa_size_t j=0; j<spectreV2.Size(); j++)
|
---|
664 | spectreV2(j) += Zmod2(cfour2(j)); // Zmod2(zp2[j]);
|
---|
665 | memcpy(procbuff+i*procpaqsz+procpaqsz/2, cfour2.Data(), sizeof(complex<r_4>)*cfour2.Size());
|
---|
666 | if (fgtimfreq) { // Remplissage tableau temps-frequence
|
---|
667 | for(sa_size_t c=1; c<timfreqV2.NCols(); c++) {
|
---|
668 | for(sa_size_t j=c*nfsmap_; j<(c+1)*nfsmap_; j++)
|
---|
669 | timfreqV2(nzm,c) += Zmod2(cfour2(j));
|
---|
670 | }
|
---|
671 | }
|
---|
672 |
|
---|
673 | // Calcul correlation (visibilite V1 * V2)
|
---|
674 | for(sa_size_t j=0; j<visiV12.Size(); j++)
|
---|
675 | visiV12(j)+=cfour1(j)*conj(cfour2(j));
|
---|
676 | // for(sa_size_t j=0; j<visiV12.Size(); j++) visiV12(j)+=zp1[j]*zp2[j];
|
---|
677 | if (nzm==0) {
|
---|
678 | spectreV1.Info()["StartFC"] = curfc;
|
---|
679 | spectreV2.Info()["StartFC"] = curfc;
|
---|
680 | visiV12.Info()["StartFC"] = curfc;
|
---|
681 | spectreV1.Info()["StartTT"] = curtt;
|
---|
682 | spectreV2.Info()["StartTT"] = curtt;
|
---|
683 | visiV12.Info()["StartTT"] = curtt;
|
---|
684 | }
|
---|
685 | nzm++;
|
---|
686 | /*----DELETE
|
---|
687 | if (nmnt==0) { xnt[0]=paq.FrameCounter(); xnt[1]=paq.TimeTag(); }
|
---|
688 | for(sa_size_t j=2700; j<2800; j++) {
|
---|
689 | ms1 += Zmod2(cfour1(j)); ms2 += Zmod2(cfour2(j));
|
---|
690 | complex<r_4> zvis = cfour1(j)*conj(cfour2(j));
|
---|
691 | ms12 += Zmod2(zvis); ms12re += zvis.real(); ms12im += zvis.imag();
|
---|
692 | ms12phi+= atan2(zvis.imag(),zvis.real());
|
---|
693 | }
|
---|
694 | nmnt++;
|
---|
695 | ----*/
|
---|
696 | totnbytesproc += paq.DataSize();
|
---|
697 | totnbytesout += (2*sizeof(complex<r_4>)*cfour1.Size());
|
---|
698 |
|
---|
699 | } // Fin de boucle sur les paquets d'une zone
|
---|
700 |
|
---|
701 | /*---- DELETE
|
---|
702 | if (nmnt>0) {
|
---|
703 | double fnorm = (double)nmnt*(2800-2700);
|
---|
704 | xnt[2] = ms1 /= fnorm;
|
---|
705 | xnt[3] = ms2 /= fnorm;
|
---|
706 | xnt[4] = ms12 /= fnorm;
|
---|
707 | xnt[5] = ms12re /= fnorm;
|
---|
708 | xnt[6] = ms12im /= fnorm;
|
---|
709 | xnt[7] = ms12phi /= fnorm;
|
---|
710 | nt.Fill(xnt);
|
---|
711 | }
|
---|
712 | ----*/
|
---|
713 | if ((nzm >= nmean_) || ((kmz==(nmax_-1))&&(nzm>1))) {
|
---|
714 | spectreV1 /= (r_4)(nzm);
|
---|
715 | spectreV2 /= (r_4)(nzm);
|
---|
716 |
|
---|
717 | // pour le calcul des moyennes et sigmas de ces spectres
|
---|
718 | moyspecV1 += spectreV1;
|
---|
719 | moyspecV2 += spectreV2;
|
---|
720 | sigspecV1 += (spectreV1 && spectreV1);
|
---|
721 | sigspecV2 += (spectreV2 && spectreV2);
|
---|
722 | nmoyspec++;
|
---|
723 |
|
---|
724 | visiV12 /= complex<r_4>((r_4)nzm, 0.);
|
---|
725 |
|
---|
726 | spectreV1.Info()["NPaqMoy"] = nzm;
|
---|
727 | spectreV2.Info()["NPaqMoy"] = nzm;
|
---|
728 | visiV12.Info()["NPaqMoy"] = nzm;
|
---|
729 | spectreV1.Info()["EndFC"] = curfc;
|
---|
730 | spectreV2.Info()["EndFC"] = curfc;
|
---|
731 | visiV12.Info()["EndFC"] = curfc;
|
---|
732 | spectreV1.Info()["EndTT"] = curtt;
|
---|
733 | spectreV2.Info()["EndTT"] = curtt;
|
---|
734 | visiV12.Info()["EndTT"] = curtt;
|
---|
735 | {
|
---|
736 | sprintf(fname,"%s_%d.ppf",path_.c_str(),(int)ifile);
|
---|
737 | POutPersist po(fname);
|
---|
738 | string tag1="specV1";
|
---|
739 | string tag2="specV2";
|
---|
740 | string tag12="visiV12";
|
---|
741 | string tagh1="tshV1";
|
---|
742 | string tagh2="tshV2";
|
---|
743 | string tagtf1="timfreqV1";
|
---|
744 | string tagtf2="timfreqV2";
|
---|
745 | if (card_==2) {
|
---|
746 | tag1 = "specV3";
|
---|
747 | tag2 = "specV4";
|
---|
748 | tagh1 = "tshV1";
|
---|
749 | tagh2 = "tshV2";
|
---|
750 | tag12="visiV34";
|
---|
751 | tagtf1="timfreqV3";
|
---|
752 | tagtf2="timfreqV4";
|
---|
753 | }
|
---|
754 | po << PPFNameTag(tag1) << spectreV1;
|
---|
755 | po << PPFNameTag(tag2) << spectreV2;
|
---|
756 | po << PPFNameTag(tag12) << visiV12;
|
---|
757 | if (fghist_) {
|
---|
758 | po << PPFNameTag(tagh1) << (*ph1);
|
---|
759 | po << PPFNameTag(tagh2) << (*ph2);
|
---|
760 |
|
---|
761 | double sspvmax[3] = {0.,0.,0.};
|
---|
762 | int_4 sspvmaxidx[3] = {-1,-1,-1};
|
---|
763 | for(int jji=1;jji<visiV12.Size()-1;jji++) {
|
---|
764 | r_4 zmv2 = Zmod2(visiV12(jji));
|
---|
765 | if (zmv2>sspvmax[2]) { sspvmax[2]=zmv2; sspvmaxidx[2]=jji; }
|
---|
766 | }
|
---|
767 | TVector<r_4>& sspv = spectreV1;
|
---|
768 | for(int ic=0; ic<2; ic++) {
|
---|
769 | if (ic==1) sspv = spectreV2;
|
---|
770 | for(int jji=1;jji<sspv.Size()-1;jji++)
|
---|
771 | if (sspv(jji)>sspvmax[ic]) { sspvmax[ic]=sspv(jji); sspvmaxidx[ic]=jji; }
|
---|
772 | if (nbsig[ic] < 1) { moysig[ic]=sigsig[ic]=-1.; }
|
---|
773 | else {
|
---|
774 | moysig[ic] /= (double)nbsig[ic];
|
---|
775 | sigsig[ic] /= (double)nbsig[ic];
|
---|
776 | sigsig[ic] -= (moysig[ic]*moysig[ic]);
|
---|
777 | sigsig[ic] = sqrt(sigsig[ic]);
|
---|
778 | cout << "===Voie " << ic << " Moy=" << moysig[ic] << " Sig=" << sigsig[ic]
|
---|
779 | << " MaxSpec Amp= " << sqrt(sspvmax[ic])/double(pq.DataSize()/2/2)
|
---|
780 | << " Pos=" << sspvmaxidx[ic] << " (NPts=" << nbsig[ic] << ")" << endl;
|
---|
781 | }
|
---|
782 | }
|
---|
783 | cout << "=== Voie1x2 MaxSpec Amp= " << sqrt(sqrt(sspvmax[2])/double(pq.DataSize()/2/2))
|
---|
784 | << " Pos=" << sspvmaxidx[2] << endl;
|
---|
785 | } // fin if (fghist_)
|
---|
786 |
|
---|
787 | if (fgtimfreq) {
|
---|
788 | timfreqV1 /= (r_4)nzm;
|
---|
789 | timfreqV2 /= (r_4)nzm;
|
---|
790 | po << PPFNameTag(tagtf1) << timfreqV1;
|
---|
791 | po << PPFNameTag(tagtf2) << timfreqV2;
|
---|
792 | }
|
---|
793 | }
|
---|
794 | spectreV1 = (r_4)(0.);
|
---|
795 | spectreV2 = (r_4)(0.);
|
---|
796 | visiV12 = complex<r_4>(0., 0.);
|
---|
797 | if (fghist_) {
|
---|
798 | ph1->Zero();
|
---|
799 | ph2->Zero();
|
---|
800 | moysig[0]=moysig[1]=0.;
|
---|
801 | sigsig[0]=sigsig[1]=0.;
|
---|
802 | nbsig[0]=nbsig[1]=0;
|
---|
803 | }
|
---|
804 | if (fgtimfreq) {
|
---|
805 | timfreqV1 = (r_4)(0.);
|
---|
806 | timfreqV2 = (r_4)(0.);
|
---|
807 | }
|
---|
808 | nzm = 0; ifile++;
|
---|
809 | // ts.SetNow();
|
---|
810 | // filog << ts << " : proc file " << fname << endl;
|
---|
811 | cout << " BRProcA2C::run() created file " << fname << card2name_(card_) << endl;
|
---|
812 | }
|
---|
813 |
|
---|
814 | memgr.FreeMemZone(mid, MemZS_ProcA);
|
---|
815 | } // Fin de boucle sur les zones a traiter
|
---|
816 | cout << " ------------ BRProcA2C::run() END " << card2name_(card_)
|
---|
817 | << " ------------ " << endl;
|
---|
818 | /*---- DELETE
|
---|
819 | {
|
---|
820 | nt.Info()["FirstTT"]=firsttt;
|
---|
821 | cout << nt;
|
---|
822 | sprintf(fname,"%s_nt.ppf",path_.c_str());
|
---|
823 | POutPersist po(fname);
|
---|
824 | po << PPFNameTag("ntv12") << nt;
|
---|
825 | cout << " BRProcA2C::run() created NTuple file " << fname << card2name_(card_) << endl;
|
---|
826 | }
|
---|
827 | ---- */
|
---|
828 | if (nmoyspec>0) { // Calcul des moyennes et sigmas des spectres
|
---|
829 | r_4 fnms = nmoyspec;
|
---|
830 | moyspecV1 /= fnms;
|
---|
831 | moyspecV2 /= fnms;
|
---|
832 | sigspecV1 /= fnms;
|
---|
833 | sigspecV2 /= fnms;
|
---|
834 | sigspecV1 -= (moyspecV1 && moyspecV1);
|
---|
835 | sigspecV2 -= (moyspecV2 && moyspecV2);
|
---|
836 | sigspecV1 = Sqrt(sigspecV1);
|
---|
837 | sigspecV2 = Sqrt(sigspecV2);
|
---|
838 | TVector<r_4> rsbV1, rsbV2; // Rapport signal/bruit
|
---|
839 | moyspecV1.DivElt(sigspecV1, rsbV1, false, true);
|
---|
840 | moyspecV2.DivElt(sigspecV2, rsbV2, false, true);
|
---|
841 | sprintf(fname,"%s_ms.ppf",path_.c_str());
|
---|
842 | POutPersist po(fname);
|
---|
843 | po << PPFNameTag("moyspecV1") << moyspecV1;
|
---|
844 | po << PPFNameTag("moyspecV2") << moyspecV2;
|
---|
845 | po << PPFNameTag("sigspecV1") << sigspecV1;
|
---|
846 | po << PPFNameTag("sigspecV2") << sigspecV2;
|
---|
847 | po << PPFNameTag("rsbV1") << rsbV1;
|
---|
848 | po << PPFNameTag("rsbV2") << rsbV2;
|
---|
849 | cout << " BRProcA2C::run() created moysigspec file " << fname << card2name_(card_) << endl;
|
---|
850 | }
|
---|
851 |
|
---|
852 | if (fghist_) {
|
---|
853 | delete ph1;
|
---|
854 | delete ph2;
|
---|
855 | }
|
---|
856 | ts.SetNow();
|
---|
857 | tm.SplitQ();
|
---|
858 | cout << " TotalProc= " << totnbytesproc/(1024*1024) << " MBytes, rate= "
|
---|
859 | << (double)(totnbytesproc)/1024./tm.PartialElapsedTimems() << " MB/s"
|
---|
860 | << " ProcDataOut=" << totnbytesout/(1024*1024) << " MB" << endl;
|
---|
861 | cout << pcheck;
|
---|
862 | cout << " BRProcA2C::run()/Timing: " << card2name_(card_) << endl;
|
---|
863 | tm.Print();
|
---|
864 | cout << " ---------------------------------------------------------- " << endl;
|
---|
865 |
|
---|
866 | }
|
---|
867 | catch (PException& exc) {
|
---|
868 | cout << " BRProcA2C::run()/catched PException " << exc.Msg() << endl;
|
---|
869 | setRC(3);
|
---|
870 | return;
|
---|
871 | }
|
---|
872 | catch(...) {
|
---|
873 | cout << " BRProcA2C::run()/catched unknown ... exception " << endl;
|
---|
874 | setRC(4);
|
---|
875 | return;
|
---|
876 | }
|
---|
877 | setRC(0);
|
---|
878 | return;
|
---|
879 | }
|
---|
880 |
|
---|
881 | //---------------------------------------------------------------------
|
---|
882 | // Classe thread de traitement 2 x 2 voies/frames (Apres BRProcA2C)
|
---|
883 | //---------------------------------------------------------------------
|
---|
884 |
|
---|
885 | /* --Methode-- */
|
---|
886 | BRProcB4C::BRProcB4C(RAcqMemZoneMgr& mem1, RAcqMemZoneMgr& mem2, string& path,
|
---|
887 | bool fgraw, uint_4 nmean, uint_4 nmax, bool fgnotrl)
|
---|
888 | : memgr1(mem1), memgr2(mem2)
|
---|
889 | {
|
---|
890 | fgraw_ = fgraw;
|
---|
891 | nmax_ = nmax;
|
---|
892 | nmean_ = nmean;
|
---|
893 | if (fgraw_) cout << " BRProcB4C::BRProcB4C() - constructeur RAW data - NMean= " << nmean_ << endl;
|
---|
894 | else cout << " BRProcB4C::BRProcB4C() - constructeur FFT data - NMean= " << nmean_ << endl;
|
---|
895 | stop_ = false;
|
---|
896 | path_ = path;
|
---|
897 | fgnotrl_ = fgnotrl;
|
---|
898 | }
|
---|
899 |
|
---|
900 | /* --Methode-- */
|
---|
901 | void BRProcB4C::Stop()
|
---|
902 | {
|
---|
903 | stop_=true;
|
---|
904 | // cout <<" BRProcB4C::Stop ... > STOP " << endl;
|
---|
905 | }
|
---|
906 |
|
---|
907 |
|
---|
908 | /* --Methode-- */
|
---|
909 | void BRProcB4C::run()
|
---|
910 | {
|
---|
911 | setRC(1);
|
---|
912 | try {
|
---|
913 | Timer tm("BRProcB4C", false);
|
---|
914 | TimeStamp ts;
|
---|
915 | BRPaqChecker pcheck1(!fgnotrl_); // Verification/comptage des paquets
|
---|
916 | BRPaqChecker pcheck2(!fgnotrl_); // Verification/comptage des paquets
|
---|
917 |
|
---|
918 | size_t totnbytesout = 0;
|
---|
919 | size_t totnbytesproc = 0;
|
---|
920 |
|
---|
921 | cout << " BRProcB4C::run() - Starting " << ts << " NMaxMemZones=" << nmax_
|
---|
922 | << " NMean=" << nmean_ << endl;
|
---|
923 | cout << " BRProcB4C::run()... - Output Data Path: " << path_ << endl;
|
---|
924 |
|
---|
925 | uint_4 paqsz = memgr1.PaqSize();
|
---|
926 | uint_4 procpaqsz = memgr1.ProcPaqSize();
|
---|
927 | if ((paqsz != memgr2.PaqSize())||(procpaqsz!= memgr2.ProcPaqSize())) {
|
---|
928 | cout << "BRProcB4C::run()/ERROR : different paquet size -> stop \n ...(PaqSz1="
|
---|
929 | << paqsz << " Sz2=" << memgr2.PaqSize() << " ProcPaqSz1="
|
---|
930 | << procpaqsz << " Sz2=" << memgr2.ProcPaqSize() << " )" << endl;
|
---|
931 | setRC(9);
|
---|
932 | return;
|
---|
933 | }
|
---|
934 |
|
---|
935 | TVector< complex<r_4> > cfour; // composant TF
|
---|
936 | BRPaquet pq(NULL, NULL, paqsz);
|
---|
937 | /*
|
---|
938 | TVector<r_4> vx(pq.DataSize()/2);
|
---|
939 | vx = (r_4)(0.);
|
---|
940 | FFTPackServer ffts;
|
---|
941 | ffts.FFTForward(vx, cfour);
|
---|
942 |
|
---|
943 | TVector< complex<r_4> > visiV13( cfour.Size() );
|
---|
944 | TVector< complex<r_4> > visiV14( cfour.Size() );
|
---|
945 | TVector< complex<r_4> > visiV23( cfour.Size() );
|
---|
946 | TVector< complex<r_4> > visiV24( cfour.Size() );
|
---|
947 | */
|
---|
948 | int szfour = pq.DataSize()/2/2+1;
|
---|
949 | // int szfour = (paqsz-40)/2+1;
|
---|
950 | TVector< complex<r_4> > visiV13( szfour );
|
---|
951 | TVector< complex<r_4> > visiV14( szfour );
|
---|
952 | TVector< complex<r_4> > visiV23( szfour );
|
---|
953 | TVector< complex<r_4> > visiV24( szfour );
|
---|
954 | // cout << " *DBG*AAAAA ---- Vectors OK" << endl;
|
---|
955 | cout << " *DBG*BRProcB4C PaqSz=" << paqsz << " ProcPaqSize=" << procpaqsz
|
---|
956 | << " procpaqsz/2=" << procpaqsz/2 << " cfour.Size()=" << szfour
|
---|
957 | << " *8=" << szfour*8 << endl;
|
---|
958 |
|
---|
959 | DataTable dt;
|
---|
960 | dt.AddLongColumn("fc1");
|
---|
961 | dt.AddLongColumn("tt1");
|
---|
962 | dt.AddLongColumn("fc2");
|
---|
963 | dt.AddLongColumn("tt2");
|
---|
964 | DataTableRow dtr = dt.EmptyRow();
|
---|
965 |
|
---|
966 | uint_4 nzm = 0;
|
---|
967 | uint_4 totnoksfc = 0;
|
---|
968 | uint_4 totnokpaq = 0;
|
---|
969 | uint_4 totnpaq = 0;
|
---|
970 | uint_4 ifile = 0;
|
---|
971 |
|
---|
972 | uint_4 curfc=0;
|
---|
973 | uint_8 curtt=0;
|
---|
974 | uint_4 curfc2=0;
|
---|
975 | uint_8 curtt2=0;
|
---|
976 | uint_8 firsttt=0;
|
---|
977 | uint_8 firsttt2=0;
|
---|
978 | bool fgfirst=true;
|
---|
979 | for (uint_4 kmz=0; kmz<nmax_; kmz++) {
|
---|
980 | uint_4 noksfc = 0;
|
---|
981 | uint_4 nokpaq = 0;
|
---|
982 | if (stop_) break;
|
---|
983 | // cout << " *DBG*BBBBB" << kmz << endl;
|
---|
984 |
|
---|
985 | int mid1 = memgr1.FindMemZoneId(MemZA_ProcB);
|
---|
986 | Byte* buff1 = memgr1.GetMemZone(mid1);
|
---|
987 | if (buff1 == NULL) {
|
---|
988 | cout << " BRProcB4C::run()/ERROR memgr.GetMemZone(" << mid1 << ") -> NULL" << endl;
|
---|
989 | break;
|
---|
990 | }
|
---|
991 | Byte* procbuff1 = memgr1.GetProcMemZone(mid1);
|
---|
992 | if (procbuff1 == NULL) {
|
---|
993 | cout << " BRProcB4C::run()/ERROR memgr.GetProcMemZone(" << mid1 << ") -> NULL" << endl;
|
---|
994 | break;
|
---|
995 | }
|
---|
996 | int mid2 = memgr2.FindMemZoneId(MemZA_ProcB);
|
---|
997 | Byte* buff2 = memgr2.GetMemZone(mid2);
|
---|
998 | if (buff1 == NULL) {
|
---|
999 | cout << " BRProcB4C::run()/ERROR memgr.GetMemZone(" << mid2 << ") -> NULL" << endl;
|
---|
1000 | break;
|
---|
1001 | }
|
---|
1002 | Byte* procbuff2 = memgr2.GetProcMemZone(mid2);
|
---|
1003 | if (procbuff2 == NULL) {
|
---|
1004 | cout << " BRProcB4C::run()/ERROR memgr.GetProcMemZone(" << mid2 << ") -> NULL" << endl;
|
---|
1005 | break;
|
---|
1006 | }
|
---|
1007 | uint_4 i1,i2;
|
---|
1008 | i1=i2=0;
|
---|
1009 | // cout << " *DBG*CCCCCC " << kmz << " memgr1.NbPaquets() =" << memgr1.NbPaquets() << endl;
|
---|
1010 | while((i1<memgr1.NbPaquets())&&(i2<memgr2.NbPaquets())) {
|
---|
1011 | BRPaquet paq1(NULL, buff1+i1*paqsz, paqsz);
|
---|
1012 | BRPaquet paq2(NULL, buff2+i2*paqsz, paqsz);
|
---|
1013 | totnpaq++;
|
---|
1014 | // cout << " DBG["<<kmz<<"] i1,i2=" << i1 <<","<<i2<<" FC1,FC2=" <<paq1.FrameCounter()
|
---|
1015 | //<<","<<paq2.FrameCounter()<<endl;
|
---|
1016 | // on ne traite que les paquets OK
|
---|
1017 | if (!pcheck1.Check(paq1)) { i1++; continue; }
|
---|
1018 | if (!pcheck2.Check(paq2)) { i2++; continue; }
|
---|
1019 | nokpaq++;
|
---|
1020 | if (paq1.FrameCounter()<paq2.FrameCounter()) { i1++; continue; }
|
---|
1021 | if (paq2.FrameCounter()<paq1.FrameCounter()) { i2++; continue; }
|
---|
1022 | // cout << " DBG["<<kmz<<"]OKOK i1,i2=" << i1 <<","<<i2<<" FC1,FC2=" <<paq1.FrameCounter()
|
---|
1023 | // <<","<<paq2.FrameCounter()<<endl;
|
---|
1024 |
|
---|
1025 | if ((i1>=memgr1.NbPaquets())||(i2>=memgr1.NbPaquets())) {
|
---|
1026 | cout << " *BUG*BUG i1=" << i1 << " i2=" << i2 << endl;
|
---|
1027 | break;
|
---|
1028 | }
|
---|
1029 | // Les deux framecounters sont identiques ...
|
---|
1030 | noksfc++;
|
---|
1031 | curfc=paq1.FrameCounter();
|
---|
1032 | curfc2=paq2.FrameCounter();
|
---|
1033 | if (fgfirst) {
|
---|
1034 | firsttt=paq1.TimeTag(); firsttt2=paq2.TimeTag();
|
---|
1035 | cout << " BRProcB4C()/Info First FC="<<curfc<<" , "<<curfc2<<" -> TT="
|
---|
1036 | << firsttt<<" , "<<firsttt2 <<endl;
|
---|
1037 | fgfirst=false;
|
---|
1038 | }
|
---|
1039 | curtt=paq1.TimeTag()-firsttt;
|
---|
1040 | curtt2=paq2.TimeTag()-firsttt2;
|
---|
1041 | dtr[0]=curfc; dtr[1]=curtt;
|
---|
1042 | dtr[2]=curfc2; dtr[3]=curtt2;
|
---|
1043 | dt.AddRow(dtr);
|
---|
1044 |
|
---|
1045 | complex<r_4>* zp1 = (complex<r_4>*)(procbuff1+i1*procpaqsz);
|
---|
1046 | complex<r_4>* zp2 = (complex<r_4>*)(procbuff1+i1*procpaqsz+procpaqsz/2);
|
---|
1047 | complex<r_4>* zp3 = (complex<r_4>*)(procbuff2+i2*procpaqsz);
|
---|
1048 | complex<r_4>* zp4 = (complex<r_4>*)(procbuff2+i2*procpaqsz+procpaqsz/2);
|
---|
1049 | for(sa_size_t j=0; j<visiV13.Size(); j++) {
|
---|
1050 | visiV13(j)+=zp1[j]*conj(zp3[j]);
|
---|
1051 | visiV14(j)+=zp1[j]*conj(zp4[j]);
|
---|
1052 | visiV23(j)+=zp2[j]*conj(zp3[j]);
|
---|
1053 | visiV24(j)+=zp2[j]*conj(zp4[j]);
|
---|
1054 | }
|
---|
1055 | if (nzm==0) {
|
---|
1056 | visiV13.Info()["StartFC"] = curfc;
|
---|
1057 | visiV14.Info()["StartFC"] = curfc;
|
---|
1058 | visiV23.Info()["StartFC"] = curfc;
|
---|
1059 | visiV24.Info()["StartFC"] = curfc;
|
---|
1060 | visiV13.Info()["StartTT"] = curtt;
|
---|
1061 | visiV14.Info()["StartTT"] = curtt;
|
---|
1062 | visiV23.Info()["StartTT"] = curtt;
|
---|
1063 | visiV24.Info()["StartTT"] = curtt;
|
---|
1064 | }
|
---|
1065 | nzm++; i1++; i2++;
|
---|
1066 | totnbytesproc += 2*paq1.DataSize();
|
---|
1067 | } // Fin de boucle sur les paquets d'une zone
|
---|
1068 | memgr1.FreeMemZone(mid1, MemZS_ProcB);
|
---|
1069 | memgr2.FreeMemZone(mid2, MemZS_ProcB);
|
---|
1070 |
|
---|
1071 | if ((nzm >= nmean_) || ((kmz==(nmax_-1))&&(nzm>1))) {
|
---|
1072 | visiV13 /= complex<r_4>((r_4)nzm, 0.);
|
---|
1073 | visiV14 /= complex<r_4>((r_4)nzm, 0.);
|
---|
1074 | visiV23 /= complex<r_4>((r_4)nzm, 0.);
|
---|
1075 | visiV24 /= complex<r_4>((r_4)nzm, 0.);
|
---|
1076 | visiV13.Info()["NPaqMoy"] = nzm;
|
---|
1077 | visiV14.Info()["NPaqMoy"] = nzm;
|
---|
1078 | visiV23.Info()["NPaqMoy"] = nzm;
|
---|
1079 | visiV24.Info()["NPaqMoy"] = nzm;
|
---|
1080 | visiV13.Info()["EndFC"] = curfc;
|
---|
1081 | visiV14.Info()["EndFC"] = curfc;
|
---|
1082 | visiV23.Info()["EndFC"] = curfc;
|
---|
1083 | visiV24.Info()["EndFC"] = curfc;
|
---|
1084 | visiV13.Info()["EndTT"] = curtt;
|
---|
1085 | visiV14.Info()["EndTT"] = curtt;
|
---|
1086 | visiV23.Info()["EndTT"] = curtt;
|
---|
1087 | visiV24.Info()["EndTT"] = curtt;
|
---|
1088 | char fname[512];
|
---|
1089 | {
|
---|
1090 | sprintf(fname,"%s_%d.ppf",path_.c_str(),(int)ifile);
|
---|
1091 | POutPersist po(fname);
|
---|
1092 | po << PPFNameTag("visiV13") << visiV13;
|
---|
1093 | po << PPFNameTag("visiV14") << visiV14;
|
---|
1094 | po << PPFNameTag("visiV23") << visiV23;
|
---|
1095 | po << PPFNameTag("visiV24") << visiV24;
|
---|
1096 | }
|
---|
1097 | visiV13 = complex<r_4>(0., 0.);
|
---|
1098 | visiV14 = complex<r_4>(0., 0.);
|
---|
1099 | visiV23 = complex<r_4>(0., 0.);
|
---|
1100 | visiV24 = complex<r_4>(0., 0.);
|
---|
1101 | nzm = 0; ifile++;
|
---|
1102 | // ts.SetNow();
|
---|
1103 | // filog << ts << " : proc file " << fname << endl;
|
---|
1104 | cout << " BRProcB4C::run() created file " << fname << endl;
|
---|
1105 | }
|
---|
1106 | double okfrac = (nokpaq>1)?((double)noksfc/(double)nokpaq*100.):0.;
|
---|
1107 | cout << "BRProcB4C ["<<kmz<<"] NOKPaq=" << nokpaq << " NSameFC=" << noksfc
|
---|
1108 | << " (" << okfrac << " %)" << endl;
|
---|
1109 | totnokpaq += nokpaq;
|
---|
1110 | totnoksfc += noksfc;
|
---|
1111 | } // Fin de boucle sur les zones a traiter
|
---|
1112 | cout << " ------------------ BRProcB4C::run() END ----------------- " << endl;
|
---|
1113 | {
|
---|
1114 | dt.Info()["FirstTT1"]=firsttt;
|
---|
1115 | dt.Info()["FirstTT2"]=firsttt2;
|
---|
1116 | cout << dt;
|
---|
1117 | char fname[512];
|
---|
1118 | sprintf(fname,"%s_fctt.ppf",path_.c_str());
|
---|
1119 | POutPersist po(fname);
|
---|
1120 | po << PPFNameTag("ttfc") << dt;
|
---|
1121 | cout << " BRProcB4C::run() created TimeTag/FrameCounter file " << fname << endl;
|
---|
1122 | }
|
---|
1123 | ts.SetNow();
|
---|
1124 | tm.SplitQ();
|
---|
1125 | cout << " TotalProc= " << totnbytesproc/(1024*1024) << " MBytes, rate= "
|
---|
1126 | << (double)(totnbytesproc)/1024./tm.PartialElapsedTimems() << " MB/s" << endl;
|
---|
1127 | double totokfrac = (totnokpaq>1)?((double)totnoksfc/(double)totnokpaq*100.):0.;
|
---|
1128 | cout << " NOkPaq1,2=" << totnokpaq << " /TotNPaq=" << totnpaq << " TotNSameFC="
|
---|
1129 | << totnoksfc << " (" << totokfrac << " %)" << endl;
|
---|
1130 | // cout << pcheck1;
|
---|
1131 | // cout << pcheck2;
|
---|
1132 | cout << " BRProcB4C::run()/Timing: \n";
|
---|
1133 | tm.Print();
|
---|
1134 | cout << " ---------------------------------------------------------- " << endl;
|
---|
1135 | }
|
---|
1136 | catch (PException& exc) {
|
---|
1137 | cout << " BRProcB4C::run()/catched PException " << exc.Msg() << endl;
|
---|
1138 | setRC(3);
|
---|
1139 | return;
|
---|
1140 | }
|
---|
1141 | catch(...) {
|
---|
1142 | cout << " BRProcB4C::run()/catched unknown ... exception " << endl;
|
---|
1143 | setRC(4);
|
---|
1144 | return;
|
---|
1145 | }
|
---|
1146 | setRC(0);
|
---|
1147 | return;
|
---|
1148 | }
|
---|
1149 |
|
---|
1150 |
|
---|