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