source: Sophya/trunk/ArchTOIPipe/ProcWSophya/simcleaner.cc@ 2351

Last change on this file since 2351 was 2054, checked in by ansari, 23 years ago

corrections pour proj ds localmap - Reza 11/6/2002
File size: 9.8 KB
Line 
1// ArchTOIPipe (C) CEA/DAPNIA/SPP IN2P3/LAL
2// Eric Aubourg
3// Christophe Magneville
4// Reza Ansari
5#include "config.h"
6
7#include "array.h"
8#include "simcleaner.h"
9#include <math.h>
10#include "toimanager.h"
11#include "pexceptions.h"
12#include "ctimer.h"
13
14#include "flagtoidef.h"
15
16SimpleCleaner::SimpleCleaner(int wsz, int nbw)
17
18{
19 SetMeanSigWindow(wsz, nbw);
20 SetCleanForMeanThrNSig();
21 SetFlaggingThrNSig();
22
23 SetRange();
24 SetFlagMask();
25
26 FillMeanSigNTuple();
27 char* noms[] = {"sn","mean","sigma"};
28 meansignt = new NTuple(3,noms);
29
30 totnscount = 0;
31 totnbblock = 0;
32 ns_outofrange = ns_flag1p = ns_flag2p = ns_flag1m = ns_flag2m = 0;
33 ns_inflagged = 0;
34 ns_cleannsig = 0;
35 gl_sum = gl_sum2 = 0.; ns_glms = 0;
36}
37
38SimpleCleaner::~SimpleCleaner()
39{
40}
41
42void SimpleCleaner::SetMeanSigWindow(int wsz, int nbw)
43{
44 mWSz = (wsz > 8) ? wsz : 8;
45 mNbW = (nbw > 3) ? nbw : 3;
46}
47
48void SimpleCleaner::PrintStatus(::ostream & os)
49{
50 os << "\n ------------------------------------------------------ \n"
51 << " SimpleCleaner::PrintStatus() - MeanWSize= " << mWSz << " mNbW="
52 << mNbW << " Range/Min=" << range_min << " Max=" << range_max << endl;
53 os << " CleanForMeanThrNSig = " << nsigclean << " NSig1= " << nsig1
54 << " NSig2 (" << min_npt2 << "pts)= " << nsig2 << endl;
55 TOIProcessor::PrintStatus(os);
56 os << " ProcessedSampleCount=" << ProcessedSampleCount()
57 << " NS_InFlagged= " << ns_inflagged
58 << " NS_OutOfRange= " << ns_outofrange
59 << " NS_CleanNSig= " << ns_cleannsig << endl
60 << " NS_Flag1(+)= " << ns_flag1p << " NS_Flag2(+)= " << ns_flag2p
61 << " NS_Flag1(-)= " << ns_flag1m << " NS_Flag2(-)= " << ns_flag2m << endl;
62 os << " NS_GlobMeanSig= " << ns_glms << " GlMean()= " << GetGlMean()
63 << " GlSigma2() " << GetGlSigma2() << endl;
64 os << " ------------------------------------------------------ " << endl;
65}
66
67void SimpleCleaner::init()
68{
69 cout << "SimpleCleaner::init" << endl;
70 declareInput("in");
71 declareOutput("out");
72 declareOutput("mean");
73 declareOutput("sigma");
74 name = "SimpleCleaner";
75}
76
77void SimpleCleaner::run()
78{
79 int snb = getMinIn();
80 int sne = getMaxIn();
81
82 bool fgout = checkOutputTOIIndex(0);
83 bool fgmean = checkOutputTOIIndex(1);
84 bool fgsig = checkOutputTOIIndex(2);
85
86 if (!checkInputTOIIndex(0)) {
87 cerr << " SimpleCleaner::run() - Input TOI (in) not connected! "
88 << endl;
89 throw ParmError("SimpleCleaner::run() Input TOI (in) not connected!");
90 }
91
92 cout << " SimpleCleaner::run() SNRange=" << snb << " - " << sne << endl;
93
94 try {
95
96 // Vecteurs pour les donnees et les sorties
97 int wsize = mWSz;
98 int hisblk = mNbW;
99 int vecsize = mWSz*mNbW;
100 Vector vinhist(vecsize);
101 TVector<uint_8> vfghist(vecsize);
102
103 Vector sumx(mNbW);
104 Vector sumx2(mNbW);
105 TVector<int_4> nbx(mNbW);
106
107 // Variables diverses
108 int k,kb,kbm,j,klast;
109 klast = snb-1;
110 totnbblock = 0;
111
112
113 int ns_flg2p_last = 0;
114 int ns_flg2m_last = 0;
115
116 // Moyenne et sigma courant
117 double cur_mean = 0.;
118 double cur_sig = 0.;
119
120 double last_sum = 0.;
121 double last_sum2 = 0.;
122
123 // Boucle sur les sampleNum
124 // 1ere partie, on traite par paquets de wsize
125 int nblk = (sne-snb+1)/wsize;
126 for(kb=0; kb<nblk; kb++) {
127 k = kb*wsize+snb;
128 kbm = kb%hisblk;
129 // Lecture d'un bloc de donnees
130 Vector vin( vinhist(Range(kbm*wsize, 0, wsize) ) );
131 TVector<uint_8> vfg( vfghist(Range(kbm*wsize, 0, wsize) ) );
132
133 getData(0, k, wsize, vin.Data(), vfg.Data());
134 double nok = 0.;
135 double sum = 0.;
136 double sum2 = 0.;
137 double cleanthrmin = range_min;
138 double cleanthrmax = range_max;
139 if ((kb > mNbW/2) && ( ns_glms > 3*wsize) ) {
140 cleanthrmin = cur_mean-nsigclean*cur_sig;
141 cleanthrmax = cur_mean+nsigclean*cur_sig;
142 }
143 for(j=0; j<wsize; j++) {
144 double x = vin(j);
145 if ((vfg(j)&flag_mask)) { ns_inflagged++; continue; }
146 if ((x > range_max) || (x < range_min) ) continue;
147 if ((x > cleanthrmax) || (x < cleanthrmin) ) {
148 ns_cleannsig++;
149 continue;
150 }
151 sum += x; sum2 += x*x; nok++;
152 }
153 if (nok > 0) {
154 gl_sum += (sum - last_sum);
155 gl_sum2 += (sum2 - last_sum2);
156 ns_glms += nok;
157 last_sum = sum; last_sum2 = sum2;
158 }
159 else {
160 sum = cur_mean; sum2 = cur_mean*cur_mean;
161 nok = 1.;
162 }
163 sumx(kbm) = sum;
164 sumx2(kbm) = sum2;
165 nbx(kbm) = nok;
166 sum = sumx.Sum();
167 sum2 = sumx2.Sum();
168 nok = nbx.Sum();
169 if (nok > 0) {
170 cur_mean = sum/(double)nok;
171 cur_sig = sum2/(double)nok-cur_mean*cur_mean;
172 cur_sig = (cur_sig > 0.) ? sqrt(cur_sig) : 0.;
173 }
174
175 if (kb < mNbW/2) continue;
176
177 Vector vinc;
178 TVector<uint_8> vfgc;
179
180 int kbs = kb-mNbW/2;
181 k = kbs*wsize+snb;
182 if (kb == nblk-1) {
183 int wszt = wsize*(kb+1-kbs);
184 vinc.ReSize(wszt);
185 vfgc.ReSize(wszt);
186 int jb = 0;
187 for(int kbsc=kbs; kbsc<nblk; kbsc++) {
188 vinc(Range(jb*wsize, 0, wsize)) =
189 vinhist(Range((kbsc%hisblk)*wsize, 0, wsize) ) ;
190 vfgc(Range(jb*wsize, 0, wsize)) =
191 vfghist(Range((kbsc%hisblk)*wsize, 0, wsize) ) ;
192 jb++;
193 }
194 wsize = wszt;
195 }
196 else {
197 vinc.Share( vinhist(Range((kbs%hisblk)*wsize, 0, wsize)) ) ;
198 vfgc.Share( vfghist(Range((kbs%hisblk)*wsize, 0, wsize)) ) ;
199 }
200
201 Vector vout(wsize);
202
203 if (fgmean) { // Sortie valeur moyenne courante
204 vout = cur_mean;
205 putData(1, k, wsize, vout.Data());
206 }
207 if (fgsig) { // Sortie sigma courante
208 vout = cur_sig;
209 putData(2, k, wsize, vout.Data());
210 }
211
212 if (fg_meansignt) { // Remplissage du NTuple mean-sigma
213 float xnt[5];
214 xnt[0] = k;
215 xnt[1] = cur_mean;
216 xnt[2] = cur_sig;
217 meansignt->Fill(xnt);
218 }
219
220 if (!fgout) {
221 klast+=wsize;
222 totnscount+=wsize;
223 totnbblock++;
224 continue;
225 }
226
227 // Calcul des flags en sortie
228 for(j=0; j<wsize; j++) {
229 double x = vinc(j);
230 if ((x > range_max) || (x < range_min) ) {
231 ns_outofrange++; vfgc(j) |= FlgToiOut;
232 }
233 if (x > cur_mean) { // Superieur a mean
234 ns_flg2m_last = 0;
235 if (x > cur_mean+nsig2*cur_sig) ns_flg2p_last++;
236 else ns_flg2p_last = 0;
237 if (x > cur_mean+nsig1*cur_sig) {
238 ns_flag1p++; vfgc(j) |= FlgToiSpike;
239 }
240 if (ns_flg2p_last == min_npt2) {
241 int jj0 = j-ns_flg2p_last+1;
242 if (jj0 < 0) jj0 = 0;
243 for(int jj = jj0; jj<=j; jj++) {
244 ns_flag2p++; vfgc(jj) |= FlgToiSpike;
245 }
246 }
247 else if (ns_flg2p_last > min_npt2) {
248 ns_flag2p++; vfgc(j) |= FlgToiSpike;
249 }
250 }
251 else { // Inferieur a mean
252 ns_flg2p_last = 0;
253 if (x < cur_mean-nsig2*cur_sig) ns_flg2m_last++;
254 else ns_flg2m_last = 0;
255 if (x < cur_mean-nsig1*cur_sig) {
256 ns_flag1m++; vfgc(j) |= FlgToiSpike;
257 }
258 if (ns_flg2m_last == min_npt2) {
259 int jj0 = j-ns_flg2m_last+1;
260 if (jj0 < 0) jj0 = 0;
261 for(int jj = jj0; jj<=j; jj++) {
262 ns_flag2m++; vfgc(jj) |= FlgToiSpike;
263 }
264 }
265 else if (ns_flg2m_last > min_npt2) {
266 ns_flag2m++; vfgc(j) |= FlgToiSpike;
267 }
268 }
269 }
270
271 putData(0, k, wsize, vinc.Data(), vfgc.Data());
272
273 klast+=wsize;
274 totnscount+=wsize;
275 totnbblock++;
276
277 } // Fin boucle sur les samples, par pas de wsize
278
279 // 3eme partie, on traite la fin du bloc d'echantillons si necessaire
280 if (klast < sne) {
281 wsize = sne-klast;
282 Vector vinc(wsize);
283 Vector vout(wsize);
284 TVector<uint_8> vfgc(wsize);
285 k = klast+1;
286 getData(0, k, wsize, vinc.Data(), vfgc.Data());
287
288 if (fgmean) { // Sortie valeur moyenne courante
289 vout = cur_mean;
290 putData(1, k, wsize, vout.Data());
291 }
292 if (fgsig) { // Sortie sigma courante
293 vout = cur_sig;
294 putData(2, k, wsize, vout.Data());
295 }
296
297 if (fg_meansignt) { // Remplissage du NTuple mean-sigma
298 float xnt[5];
299 xnt[0] = k;
300 xnt[1] = cur_mean;
301 xnt[2] = cur_sig;
302 meansignt->Fill(xnt);
303 }
304
305 if (fgout) for(j=0; j<wsize; j++) {
306 double x = vinc(j);
307 if ((vfgc(j)&flag_mask)) ns_inflagged++;
308 if ((x > range_max) || (x < range_min) ) {
309 ns_outofrange++; vfgc(j) |= FlgToiOut;
310 }
311 if (x > cur_mean) { // Superieur a mean
312 ns_flg2m_last = 0;
313 if (x > cur_mean+nsig2*cur_sig) ns_flg2p_last++;
314 else ns_flg2p_last = 0;
315 if (x > cur_mean+nsig1*cur_sig) {
316 ns_flag1p++; vfgc(j) |= FlgToiSpike;
317 }
318 if (ns_flg2p_last == min_npt2) {
319 int jj0 = j-ns_flg2p_last+1;
320 if (jj0 < 0) jj0 = 0;
321 for(int jj = jj0; jj<=j; jj++) {
322 ns_flag2p++; vfgc(jj) |= FlgToiSpike;
323 }
324 }
325 else if (ns_flg2p_last > min_npt2) {
326 ns_flag2p++; vfgc(j) |= FlgToiSpike;
327 }
328 }
329 else { // Inferieur a mean
330 ns_flg2p_last = 0;
331 if (x < cur_mean-nsig2*cur_sig) ns_flg2m_last++;
332 else ns_flg2m_last = 0;
333 if (x < cur_mean-nsig1*cur_sig) {
334 ns_flag1m++; vfgc(j) |= FlgToiSpike;
335 }
336 if (ns_flg2m_last == min_npt2) {
337 int jj0 = j-ns_flg2m_last+1;
338 if (jj0 < 0) jj0 = 0;
339 for(int jj = jj0; jj<=j; jj++) {
340 ns_flag2m++; vfgc(jj) |= FlgToiSpike;
341 }
342 }
343 else if (ns_flg2m_last > min_npt2) {
344 ns_flag2m++; vfgc(j) |= FlgToiSpike;
345 }
346 }
347 putData(0, k, wsize, vinc.Data(), vfgc.Data());
348 }
349
350
351 klast+=wsize;
352 totnscount+=wsize;
353 totnbblock++;
354 }
355
356 meansignt->Info()["SampleCount"] = ProcessedSampleCount();
357 meansignt->Info()["GlMean"] = GetGlMean();
358 meansignt->Info()["GlSigma2"] = GetGlSigma2();
359 meansignt->Info()["OutOfRange"] = ns_outofrange;
360 meansignt->Info()["NSFlag1+"] = ns_flag1p;
361 meansignt->Info()["NSFlag1-"] = ns_flag1m;
362 meansignt->Info()["NSFlag2+"] = ns_flag2p;
363 meansignt->Info()["NSFlag2-"] = ns_flag2m;
364
365 cout << " SimpleCleaner::run() - End of processing "
366 << " TotNbBlocks= " << totnbblock << " ProcSamples=" << totnscount << endl;
367
368 } // Bloc try
369
370 catch (PException & exc) {
371 cerr << "SimpleCleaner::run() Catched Exception " << (string)typeid(exc).name()
372 << "\n .... Msg= " << exc.Msg() << endl;
373 }
374
375}
Note: See TracBrowser for help on using the repository browser.