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

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

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