source: Sophya/trunk/ArchTOIPipe/ProcWSophya/simoffset.cc@ 2000

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

Corrections diverses et introduction d'une classe de calcul de ligne de
base par ajustement glissant d'un polynome (SimpleOffsetEstimator)
avec son programme test - Reza 14/5/2002
File size: 7.2 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 "simoffset.h"
9#include <math.h>
10#include "toimanager.h"
11#include "pexceptions.h"
12#include "ctimer.h"
13#include "xntuple.h"
14
15#include "flagtoidef.h"
16
17SimpleOffsetEstimator::SimpleOffsetEstimator(int mwsz, int nptfit, int degpol)
18 : poly(degpol)
19{
20 mWSz = (mwsz > 8) ? mwsz : 8;
21 nPtFit = (nptfit > degpol+2) ? nptfit : degpol+2;
22 totnscount = 0;
23 totnbblock = 0;
24 SavePolyNTuple();
25}
26
27SimpleOffsetEstimator::~SimpleOffsetEstimator()
28{
29}
30
31void SimpleOffsetEstimator::PrintStatus(::ostream & os)
32{
33 os << "\n ------------------------------------------------------ \n"
34 << " SimpleDeglitcher::PrintStatus() - MeanWSize= " << mWSz << " NPtFit="
35 << nPtFit << " DegPoly=" << poly.Degre() << endl;
36 TOIProcessor::PrintStatus(os);
37 os << " ProcessedSampleCount=" << ProcessedSampleCount() << endl;
38 os << " ------------------------------------------------------ " << endl;
39}
40
41void SimpleOffsetEstimator::init()
42{
43 cout << "SimpleOffsetEstimator::init" << endl;
44 declareInput("in");
45 declareOutput("offset");
46 declareOutput("out");
47 declareOutput("incopie");
48 declareOutput("poly_a0");
49 declareOutput("poly_a1");
50 declareOutput("poly_a2");
51 declareOutput("poly_sn0");
52 declareOutput("mean_y");
53 declareOutput("sig_y");
54 declareOutput("mean_x");
55 name = "SimpleOffsetEstimator";
56}
57
58void SimpleOffsetEstimator::run()
59{
60 int snb = getMinIn();
61 int sne = getMaxIn();
62
63 bool fgoffset = checkOutputTOIIndex(0);
64 bool fgout = checkOutputTOIIndex(1);
65 bool fgincopie = checkOutputTOIIndex(2);
66 bool fga0 = checkOutputTOIIndex(3);
67 bool fga1 = checkOutputTOIIndex(4);
68 bool fga2 = checkOutputTOIIndex(5);
69 bool fgsn0 = checkOutputTOIIndex(6);
70 bool fgmeany = checkOutputTOIIndex(7);
71 bool fgsigy = checkOutputTOIIndex(8);
72 bool fgmeanx = checkOutputTOIIndex(9);
73
74 if (!checkInputTOIIndex(0)) {
75 cerr << " SimpleOffsetEstimator::run() - Input TOI (in) not connected! "
76 << endl;
77 throw ParmError("SimpleOffsetEstimator::run() Input TOI (in) not connected!");
78 }
79 if (!fgoffset && !fgout) {
80 cerr << " SimpleOffsetEstimator::run() - No Output TOI (offset/in-offset) connected! "
81 << endl;
82 throw ParmError(" SimpleOffsetEstimator::run() No output TOI (offset/in-offset) connected!");
83 }
84
85 cout << " SimpleOffsetEstimator::run() SNRange=" << snb << " - " << sne << endl;
86
87 // NTuple pour sauvegarde des coeff de poly
88 char * nomsnt[] = {"sncur", "sn0", "meanx", "meany", "sigy", "a0", "a1", "a2", "ycur"};
89 XNTuple xntp(0, 9, 0, 0, nomsnt);
90
91 try {
92
93 // Vecteurs pour les donnees et les sorties
94 int wsize = mWSz;
95
96 Vector vin(wsize);
97 Vector voff(wsize);
98 Vector vout(wsize);
99 TVector<uint_8> vfg(wsize);
100
101 // Pour le fit
102 Vector errCoef(3);
103
104 Vector X(nPtFit);
105 Vector X0(nPtFit);
106 Vector Y(nPtFit);
107 Vector YErr(nPtFit);
108
109 // Variables diverses
110 int k,i,j,klast;
111 int nks = 0;
112 klast = snb-1;
113 totnbblock = 0;
114
115
116 int nbblkok = 0;
117
118 double sn0 = 0.;
119 double nok = 0.;
120 double mean = 0.;
121 double sig = 0.;
122 double meanx = 0.;
123
124 // Boucle sur les sampleNum
125 // 1er partie, on traite par paquets de wsize
126
127 for(k=snb;k<=sne-wsize+1;k+=wsize) {
128 // Lecture d'un bloc de donnees
129 getData(0, k, wsize, vin.Data(), vfg.Data());
130
131 // Calcul moyenne et sigma du bloc
132 nok = 0.; meanx = 0.;
133 mean = 0.; sig = 0.;
134
135 for(j=0; j<wsize; j++) {
136 if ( vfg(j) ) continue;
137 mean += vin(j);
138 sig += vin(j)*vin(j);
139 meanx += k+j;
140 nok++;
141 }
142 if (nbblkok == 0) {
143 X = RegularSequence(k+wsize*0.5, (double)wsize);
144 Y = mean;
145 YErr = (nok > 0.5) ? sqrt(mean) : 1.;
146 }
147
148 sn0 = (double)(k+wsize/2);
149
150 if (nok > 3.) {
151 mean /= nok;
152 meanx /= nok;
153 sig = sig/nok-mean*mean;
154 int kk = nbblkok%nPtFit;
155 nbblkok++;
156 Y(kk) = mean;
157 YErr(kk) = sig;
158 X(kk) = meanx;
159 }
160
161 X0 = X;
162 X0 -= sn0;
163 poly.Fit(X,Y,YErr,poly.Degre(),errCoef);
164
165
166 // Calcul des valeurs d'offset en sortie
167 for(j=0; j<wsize; j++)
168 voff(j) = poly(k+j-sn0);
169
170 if (fgoffset) putData(0, k, wsize, voff.Data());
171 if (fgincopie) putData(2, k, wsize, vin.Data(), vfg.Data());
172 if (fgout) {
173 vin -= voff;
174 putData(1, k, wsize, vin.Data(), vfg.Data());
175 }
176
177 if (fga0) {
178 vout = poly[0];
179 putData(3, k, wsize, vout.Data());
180 }
181 if (fga1) {
182 vout = poly[1];
183 putData(4, k, wsize, vout.Data());
184 }
185 if (fga2) {
186 vout = poly[2];
187 putData(5, k, wsize, vout.Data());
188 }
189 if (fgsn0) {
190 vout = sn0;
191 putData(6, k, wsize, vout.Data());
192 }
193 if (fgmeany) {
194 vout = mean;
195 putData(7, k, wsize, vout.Data());
196 }
197 if (fgsigy) {
198 vout = sig;
199 putData(8, k, wsize, vout.Data());
200 }
201
202 if (fgmeanx) {
203 vout = meanx;
204 putData(9, k, wsize, vout.Data());
205 }
206
207 if (ntpoly) { // Remplissage du XNTuple de controle
208 char * nomsnt[] = {"sncur", "sn0", "meanx", "meany", "sigy", "a0", "a1", "a2", "ycur"};
209 float xnt[10];
210 xnt[0] = k;
211 xnt[1] = sn0;
212 xnt[2] = meanx;
213 xnt[3] = mean;
214 xnt[4] = sig;
215 xnt[5] = poly[0];
216 xnt[6] = poly[1];
217 xnt[7] = poly[2];
218 xnt[8] = poly(k-sn0);
219 xntp.Fill(NULL, xnt, NULL, NULL);
220 }
221
222 klast+=wsize;
223 totnscount+=wsize;
224 totnbblock++;
225
226 } // Fin boucle sur les samples, par pas de wsize
227
228 // 2eme partie, on traite la fin du bloc d'echantillons si necessaire
229 if (klast < sne) {
230 wsize = sne-klast;
231 vin.ReSize(wsize);
232 voff.ReSize(wsize);
233 vout.ReSize(wsize);
234 vfg.ReSize(wsize);
235 getData(0, k, wsize, vin.Data(), vfg.Data());
236 for(j=0; j<wsize; j++)
237 voff(j) = poly(k+j-sn0);
238 if (fgoffset) putData(0, k, wsize, voff.Data());
239 if (fgincopie) putData(2, k, wsize, vin.Data(), vfg.Data());
240 if (fgout) {
241 vin -= voff;
242 putData(1, k, wsize, vin.Data(), vfg.Data());
243 }
244
245 if (fga0) {
246 vout = poly[0];
247 putData(3, k, wsize, vout.Data());
248 }
249 if (fga1) {
250 vout = poly[1];
251 putData(4, k, wsize, vout.Data());
252 }
253 if (fga2) {
254 vout = poly[2];
255 putData(5, k, wsize, vout.Data());
256 }
257 if (fgsn0) {
258 vout = sn0;
259 putData(6, k, wsize, vout.Data());
260 }
261 if (fgmeany) {
262 vout = mean;
263 putData(7, k, wsize, vout.Data());
264 }
265 if (fgsigy) {
266 vout = sig;
267 putData(8, k, wsize, vout.Data());
268 }
269
270 if (fgmeanx) {
271 vout = meanx;
272 putData(9, k, wsize, vout.Data());
273 }
274
275 klast+=wsize;
276 totnscount+=wsize;
277 totnbblock++;
278 }
279
280 cout << " SimpleOffsetEstimator::run() - End of processing "
281 << " TotNbBlocks= " << totnbblock << " ProcSamples=" << totnscount << endl;
282
283 } // Bloc try
284
285 catch (PException & exc) {
286 cerr << "SimpleOffsetEstimator::run() Catched Exception " << (string)typeid(exc).name()
287 << "\n .... Msg= " << exc.Msg() << endl;
288 }
289
290 if (ntpoly) {
291 if (ntpolyname.length() < 1) ntpolyname = "simoffset.ppf";
292 POutPersist pos(ntpolyname);
293 cout << " SimpleOffsetEstimator::run()/Info : Writing poly ntuple to PPF file " << ntpolyname << endl;
294 pos << xntp;
295 }
296}
Note: See TracBrowser for help on using the repository browser.